Shack-Hartmann Wavefront Sensor Laser Tomography

In [1]:
import numpy as np
import ceo
%pylab inline

Populating the interactive namespace from numpy and matplotlib

A $N_L\times N_L$ Shack-Hartmann wavefront sensor (WFS) on a $D$ diameter telescope is going to be modeled with $n_P\times n_P$ pixels per lenslet on the detector.

In [2]:
nLenslet = 60
NA = nLenslet + 1;
D = 25.5
n = 6
nPx = n*nLenslet + 1

A arbitrary number of guide stars (GS) can be used by one WFS model, the assumption beeing that all the GS uses exactly the same WFS. Here the GS are randomly distributed in a 2arcmin diameter field of regard.

In [3]:
nGS = 6
gs = ceo.Source("K",
                 zenith=np.ones(nGS)*30*math.pi/180/3600,
                 azimuth=np.linspace(0,nGS-1,nGS)*2*math.pi/nGS,
                 height = 90e3,
                 resolution=(nPx,nPx))
calib_src = ceo.Source("K", resolution=(nPx,nPx))

Next the telescope pupil mask is defined and applied to the GSs.

In [4]:
tel = ceo.GMT(nPx,D)
gs.masked(tel)
calib_src.masked(tel)

The shack-Hartmann WFS is defined with the shackHartmann class. The mandatory parameters of the class constructor are:

  • the lenslet array size $N_L$,
  • the number of pixel per lenslet in the pupil plane $n_P$ with pixels on the lenslet edges, meaning that the wavefront of the corresponding GSs must be sampled with $(N_L n_P+1) \times (N_L n_P+1)$ pixels,

  • the lenslet pitch $d$.

    By default, the WFS imagelets are Nyquist sampled meaning that the default value of the discrete Fourier transform over-sampling factor is set to 2, i.e. $\alpha\equiv$ DFT_osf=2. The size of each imagelet is then $\alpha n_P \times \alpha n_P$.

    The detector framelet sizes are N_PX_IMAGE$\times $ N_PX_IMAGE pixels each with the default value N_PX_IMAGE $=n_P$. If the imagelets are larger, they are cropped to the specified size, if they are smaller, the framelets are padded with zeros.

    Finally the framelets can be binned down by a factor $\beta\equiv$ BIN_IMAGE leading to a final framelet size of $${N_I\over\beta}\times{N_I\over\beta}$$ with $N_I\equiv$ N_PX_IMAGE. The default value of BIN_IMAGE is 1.

    The pixel scale is given by: $${\beta\lambda \over \alpha d}$$ and the lenslet field-of-view is $$N_I{\lambda\over \alpha d}$$ with $\lambda$ the wavelength.

    If more than on GS is assigned to the WFS, the parameter N_GS needs to be set.

In [5]:
d = D/nLenslet
wfs = ceo.ShackHartmann(nLenslet, n, d, N_PX_IMAGE=2*(n+1))
wfs.calibrate(calib_src,0.5)
imshow(wfs.flux.host().T,interpolation='none')
Out[5]:
<matplotlib.image.AxesImage at 0x2b368c38ad50>
In [6]:
px_scale = 2.179e-6/d/2
coef_med = []
px = arange(0,1,0.05)
for k in px:
    wfs.pointing(np.array(-px_scale*k,dtype=np.float32,ndmin=1),np.array(0.0,dtype=np.float32,ndmin=1))
    wfs.reset()
    wfs.analyze(calib_src)
    c = wfs.c.host()
    cx = c[0,0:c.size/2]
    m = wfs.valid_lenslet.f.host()
    cx = cx[m.flatten()>0]
    #print k, np.mean(cx/px_scale/k) , np.median(cx/px_scale/k)
    coef_med.append(np.median(cx/px_scale))
cp = np.polyfit(px, coef_med, 1)
print cp
plot(px,coef_med,px,coef_med/cp[0])
grid()

[ 1.00597442 -0.02863912]
In [7]:
wfs = ceo.ShackHartmann(nLenslet, n, d, N_GS = nGS, N_PX_IMAGE=2*(n+1))
#wfs.slopes_gain = 1.0/cp[0]

The WFS reference slopes and valid lenslets are set with the calibrate method passing a Source object which wavefront sets the reference slopes and the lenslet intensity threshold used to discard the lenset with too litle illumination.

In [8]:
wfs.calibrate(gs,0.5)
In [9]:
print wfs.valid_actuator.f.host().shape

(22326, 1)
In [10]:
validActuator = wfs.valid_actuator
validActuator_f = validActuator.f
imshow(validActuator_f.host(shape=((nLenslet+1)*nGS,(nLenslet+1))).T,interpolation='None')
validActuator_f.host_data.sum()
Out[10]:
13740.0

The source wavefront can now be analyzed with respect to the former calibration.

In [11]:
wfs.analyze(gs)

The detector frame is saved in the cuFloatArray object frame.

In [12]:
#figure(figsize=(12,12))
#imshow(wfs.frame.host().transpose(),interpolation='none')

A Atmosphere object is defined next and the GSs are propagated through:

In [13]:
atm =ceo.GmtAtmosphere(20e-2,30)
p = D/(nPx-1)
atm.get_phase_screen(gs,  p, nPx, p, nPx, 0.0)
In [14]:
#figure(figsize=(12,12))
#imshow(gs.phase.host().transpose(),interpolation='none')

The WFS detector need to be reset before proceeding with a new wavefront analysis.

In [15]:
wfs.reset()
wfs.analyze(gs)
In [16]:
#figure(figsize=(12,12))
#imshow(wfs.frame.host().transpose(),interpolation='none')

The WFS centroids are accessible through the cuFloatArray object c.

In [17]:
#figure(figsize=(12,12))
#imshow(wfs.c.host().reshape(2*nLenslet*nGS
#                            ,nLenslet).transpose(),interpolation='none')
#colorbar()
In [18]:
src = ceo.Source("K",resolution=(NA,NA))
src.masked(validActuator)
In [19]:
lmmse = ceo.LmmseSH(atm,gs,src,wfs,"MINRES")
In [20]:
lmmse.estimation(wfs)
In [21]:
lmmse_phase = lmmse.phase
mask_actuator = validActuator_f.host_data[0:nLenslet+1,:]
imshow(lmmse_phase.host(units='micron',zm=True,mask=mask_actuator),interpolation='none')
colorbar()
Out[21]:
<matplotlib.colorbar.Colorbar at 0x2b368c695610>
In [22]:
atm.get_phase_screen(src,d,NA,d,NA,0.0)
In [23]:
src_phase = src.phase
ps_e = src_phase.host(units='micron',zm=True,mask=mask_actuator) - \
    lmmse_phase.host(units='micron',zm=True,mask=mask_actuator)
wfe_rms = np.std(ps_e[mask_actuator!=0])*1e3
print "wavefront error: %6.2fnm" % wfe_rms
imshow(np.concatenate((src_phase.host_data,lmmse_phase.host_data),axis=1),interpolation='none')
colorbar()

wavefront error: 134.45nm
Out[23]:
<matplotlib.colorbar.Colorbar at 0x2b368c803e10>
In [24]:
imshow(ps_e*1e3,interpolation='None')
colorbar()
Out[24]:
<matplotlib.colorbar.Colorbar at 0x2b368c967b90>
In [25]:
np.exp(-(2*math.pi*wfe_rms/2.2e3)**2)
Out[25]:
0.86291018127947128
In [26]: