Implement KolmogorovWFE

[josePhoenix]

See

http://www.opticsinfobase.org/view_article.cfm?gotourl=http%3A%2F%2Fwww%2Eopticsinfobase%2Eorg%2FDirectPDFAccess%2F8E2A4176%2DED0A%2D7994%2DFB0AC49CECB235DF%5F142887%2Epdf%3Fda%3D1%26id%3D142887%26seq%3D0%26mobile%3Dno&org=

http://optics.nuigalway.ie/people/chris/chrispapers/Paper066.pdf

(Captured from docstring of placeholder class removed from wfe.py)

[mperrin]

Here's some code I wrote last year in the course of gpipsfs development. Basic Kolmogorov phase screen, though limited in its ability to properly model turbulence at larger angular scales (outer working angle).

This would need work to turn into an OpticalElement class. For now I'm just pasting it here so I don't forget about it whenever someone eventually gets around to working on this.

# Based on IDL code in fgui.pro by Tuan Do
r0=0.18 # meters
shape=(200,200)
pixelscale=0.05 #meters/pixel

# set up indices arrays
dk = 1./(np.asarray(shape)*pixelscale)
y,x = np.indices(shape)
y=(y-shape[0]/2.)*dk[0]
x=(x-shape[1]/2.)*dk[1]
ksq=x**2+y**2
ksq[shape[0]/2,shape[1]/2]=1

#kolmogorov PSD
psd=0.023*(2*np.pi)**(5./6) * r0**(-5./6) * ksq**(-11/6) * np.sqrt(dk[0]*dk[1])
psd[shape[0]/2,shape[1]/2]=0  # set piston to 0

# Generate complex independent, Gaussian, random numbers with zero mean and unit variance
w_r=np.random.normal(size=shape)
w_i=np.random.normal(size=shape)
w=w_r + np.complex(0,1)*w_i

#multiply by sqrt of PSD
wf=w*np.sqrt(psd)
out = np.fft.fft2(np.fft.fftshift(wf)).real

And here's some code for demonstrating that the output follows the 5/3 power law, but only up to some finite spatial scale:

def structure_fn(array,npts=150):
    """ Two point correlation function
    quick and hacky calculaion just using cardinal directions
    """
    ans=np.zeros((npts+1))
    for i in range(1,npts+1):
        for shift,axis in ((1,0),(-1,0),(1,1),(-1,1)):
            ans[i]=((array-np.roll(array,shift*i,axis=axis))**2).mean()
        ans[i]/4
    return ans
res = structure_fn(out)
plt.loglog(res)
x=np.arange(40)
plt.plot(x,res[1]*(x/x[1])**(5./3),ls=":")
print res[0]*(x)**(5./3)

[mperrin]

This code is now available in poppy-ized form in the kolmogorov branch. warning, not yet ready for prime time. Something's scaled incorrectly somewhere so I'm getting WFE in excess of 1 meter in the output arrays. I don't have time to debug this now but wanted to make the code available in case others have time to work on it.

[mperrin]

See #249 which is a much better version of this.