Do-loop conversion

*To*: mathgroup at smc.vnet.net*Subject*: [mg50211] Do-loop conversion*From*: "Rich Matyi" <rjmatyi at comcast.net>*Date*: Sat, 21 Aug 2004 03:04:12 -0400 (EDT)*Sender*: owner-wri-mathgroup at wolfram.com

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: Do[ 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 + (fc*prior*pf^2)); 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 comcast.net

**Follow-Ups**:**Re: Do-loop conversion***From:*DrBob <drbob@bigfoot.com>

**Re: Re: Re: FindMinimum and the minimum-radius circle**

**Re: Re: Plotting a contour plot with cylindrical co-ordinates**

**PDE with boundary condition ODE**

**Re: Do-loop conversion**