Laser Tomography

Phase screens are computed with an Atmosphere and a Source object. Lets import the ceo module first.

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

Populating the interactive namespace from numpy and matplotlib

The atmosphere is defined with

In [2]:
atm =  ceo.GmtAtmosphere(0.15,60)

A constellation of 6 Laser guide stars evenly located on a 1 arcmin diameter circle is defined first,

In [3]:
NL = 60
NA = NL+1
lgs = ceo.Source("V",
                 zenith=np.ones(6)*30*math.pi/180/3600,
                 azimuth=np.linspace(0,5,6)*2*math.pi/6,
                 height = 90e3,
                 resolution=(NA,NA))

The telescope pupil is defined as the Giant Magellan Telescope and mask for the deformable mirror actuators is also set.

In [4]:
D = 25.5
#tel = ceo.Telescope(NL*16)
#dm  = ceo.Telescope(NA)
tel = ceo.GMT(NL*16,D)
dm  = ceo.Mask(NA,D)

A Centroiding object is defined, it will contains the phase screen gradient corresponding to the LGS. The fried_geometry method computes the DM valid actuator mask according to the telescope pupil shape and the given intensity threshlod. The lgs Source object is masked with the dm mask.

In [5]:
d = D/NL
cog = ceo.Centroiding(NL,N_SOURCE=lgs.size)
cog.fried_geometry(dm, tel, 16, 0.5)
lgs.masked(dm)
In [6]:
dm_mask = dm.f

The phase screen gradient is computed with the Atmosphere method get_phase_screen_gradient. The gradient is computed over a square lenslet array of size $N_L \times N_L$ with $d$ the pitch in meter. The phase screen gradient is computed for a given Source object that contains one or more guide stars. The phase screen gradient can be computed for a given time delay. The $c_x$ and $c_y$ centroids are saved in a Centroiding object.

In [7]:
atm.get_phase_screen_gradient(cog,NL,d,lgs,0.0)
c = cog.c.host(units='arcsec')
In [8]:
imshow(c.reshape(NL*lgs.size*2,NL).transpose(),interpolation='none')
#ceog.heatmap(c.reshape(NL*6*2,NL).transpose(), filename=PLOTLY_PATH+"wavefront gradient")
Out[8]:
<matplotlib.image.AxesImage at 0x2b9737833a10>

The on-axis source is defined and is propagated throught the atmosphere:

In [9]:
src = ceo.Source("K",resolution=(NA,NA))
src.masked(dm)
atm.get_phase_screen(src,d,NA,d,NA,0.0)

From the phase gradient, the phase screen can be reconstructed with a linear minimim mean square error reconstructor (LMMSE). A Lmmse object is used to perform the phase estimation. The parameters are:

  • an Atmosphere object,
  • a Source object representing the guide star(s),
  • a Source object representing the star(s) in the estimation direction(s),
  • the wavefront sampling step in meter,
  • the number of sample across the wavefront,
  • a Mask object representing the pupil,
  • the iterative solver.
In [10]:
et = ceo.StopWatch()
et.tic()
src_lmmse = ceo.Lmmse(atm,lgs,src,d,NL,dm,"MINRES")
et.toc()
print "ET = %.2fms"%et.elapsedTime

ET = 23625.15ms
In [11]:
et.tic()
src_lmmse.reset()
src_lmmse.estimation(cog)
et.toc()
src_phase = src.phase
src_lmmse_phase = src_lmmse.phase
print "ET = %.2fms"%et.elapsedTime

ET = 59.54ms
In [12]:
ps_e = src_lmmse_phase.host(units='micron',
                            zm=True,mask=dm_mask.host()) - src_phase.host(units='micron',zm=True,mask=dm_mask.host_data)
print "wavefront error: %6.2fnm" % (np.std(ps_e[dm_mask.host_data.reshape((NA,NA))!=0])*1e3)

wavefront error: 129.02nm
In [13]:
imshow(np.concatenate((src_phase.host_data, src_lmmse_phase.host_data),axis=1),
             interpolation='none')
colorbar()
Out[13]:
<matplotlib.colorbar.Colorbar at 0x2b97379b5190>
In [14]:
imshow(ps_e*1e3,interpolation='none')
colorbar()
Out[14]:
<matplotlib.colorbar.Colorbar at 0x2b9737b21990>
In [15]: