Reflection spectrum from a Fibonacci quasicrystal

Warning: this question may come across as too specific. Maybe it would help if I set values for the probability parameters below - they'd be based on the Fresnel equations. But anyway, I think this is an interesting problem - or at least, I'd like to know how to find the answer.

The Fibonacci word is formed by concatenation in the same manner that the Fibonacci sequence is formed by addition, thus:

W1 = 0, W2 = 01, W3 = W2W1 = 010, Wn = Wn-1Wn-2

W = 0100101001001010010100100101001001...

A quasicrystal is a quasiperiodic structure. One of the simplest is the Fibonacci quasicrystal, a one-dimensional structure which can be defined by replacing each 0 in W with one type of optically transparent layer and each 1 with another type. Partial reflection occurs at the interface between each pair of adjacent layers, with a probability dependent on the types of the layers. Light which is not reflected passes instead into the next layer.

So we have five parameters determining the probabilities of reflection when passing from one medium into another: p (air into 0), q (0 into 1), r (1 into 0), s (0 into 0), t (0 into air).

Light of wavelength λ undergoes a phase change of Λ/λ after passing through a 1 layer, and φΛ/λ after passing through a 0 layer, where φ is the golden ratio (1+√5)/2.

Finally, light reflected from the quasicrystal interferes with itself so that the resulting waveform for wavelength λ is

f(τ) = p cos(2πcτ/λ) + (1-p)q(1-t) cos(2(πcτ+φΛ)/λ) + (1-p)(1-q)r(1-r)(1-t) cos(2(πcτ+(φ+1)Λ)/λ) + ...

and the associated intensity is the root mean square of f over the entire time domain.

Pretty horrific, I know. But (in theory at least) it can all be written in terms of trig functions of πcτ, φΛ/λ and Λ/λ along with sums of products of the probabilities. The question is, how can I actually evaluate the intensity associated to λ? Is it even vaguely possible without lots of numerical computation?

-
It looks like $t$ has two different meaning above? – Stopple Jul 19 '11 at 19:46
Sorry, you're quite right. Relabelled one of them to τ accordingly. – Robin Saunders Jul 19 '11 at 20:54
Did you try plotting the partial sums of f(τ) to see a bit what it looks like? – André Henriques Jul 19 '11 at 20:59
First I need a good way to work out the coefficients of each of the waveform components. To simplify, let's suppose p = q = r = s = t. We can carry out the following loop: select a value c from Geom(p) (the sign of c should alternate for each pass of the loop); add c to the current position k; find the current depth, ϕ*floor[(k+1)/ϕ] + floor[(k+1)/ϕ²]; find the unsigned change in depth; add to total distance travelled. Stop the loop when the position becomes negative. The result is a random variable supported on {mϕ+n : m, n natural numbers}. But what is its probability distribution? – Robin Saunders Jul 25 '11 at 14:44