""" Module for various Adaptive Optics (AO) simulations.
A key reference (among many) is `Guyon (2005) <http://adsabs.harvard.edu/abs/2005ApJ...629..592G>`_
To compute Guyon's C0-C6 terms, use the following functions:
C0, C1 -- :func:`contrast_uncorrected`
C2, C3 -- :func:`contrast_residual`
C4, C5, C6 -- :func:`contrast_chromaticity`
A quick test is just to run :func:`run_example`; look under the hood
at its source code to see more details.
:DEPENDENCIES:
:func:`nsdata.nAir` -- refractive index of air (for :func:`contrast_chromaticity`)
PyLab -- for plotting the figure in :func:`run_example`
:TO DO:
Add functions for beta_p (wavefront sensor sensitivity) for various
sensor technologies: Shack-Hartmann, Curvature, etc.
:NOTE:
Note that the coefficient "0.484" in Eqs. 16 and 17 of Guyon 2005
has been corrected here to its intended value of 0.22*0.22=0.0484
(O. Guyon, private communication, Jan 2013).
"""
#import pdb
import numpy as np
from nsdata import nAir
[docs]def run_example():
"""Run an example simulation, and reproduce Fig. 12 of Guyon 2005.
See the source code of this function for details.
Note that as published, the coefficients in Guyon 2005 Eqs. 16-17
are incorrect (they should be 0.0484, not 0.484). This function
uses the correct coefficient, and so differs from the published
version of Fig. 12.
"""
# 2012-12-05 21:08 IJMC: Created
import pylab as py
D = 8
lam_sci = 1.6e-6
lam_wfs = 0.55e-6
F = 9.74e7 # Because m_V = 5
beta_p = 1.0
beta_a = 1.0
r0 = 0.2
lam_r0 = 0.5e-6
v = 10
z, dz, cn2 = makecn2()
alpha_arcsec = np.logspace(-1, 1, 200) # angular separation, in arcsec
X, Y, dX, dY = calc_XY(z, dz, cn2, alpha_arcsec, lam_sci, lam_wfs=lam_wfs)
C0, C1 = contrast_uncorrected(lam_r0, r0, z, dz, cn2, alpha_arcsec, \
lam_sci, D, X=X, Y=Y)
C2, C3 = contrast_residual(lam_r0, r0, z, dz, cn2, v, alpha_arcsec, F, \
beta_p, beta_a, lam_wfs, lam_sci, D, \
X=X, Y=Y)
C4, C5, C6 = contrast_chromaticity(lam_r0, r0, z, dz, cn2, alpha_arcsec, \
lam_wfs, lam_sci, D, \
X=X, Y=Y, dX=dX, dY=dY, C0=C0, C1=C1)
py.figure()
contrasts = [C0, C1, C2, C3, C4, C5, C6]
[py.loglog(alpha_arcsec, C, linewidth=2) for C in contrasts]
py.legend(['C%i' % ii for ii in range(7)]+['Phase only', 'Phase+Amplitude'], 1)
py.axis([.1,10,1e-11,1e-3])
py.xlabel('Angular Separation [arcsec]')
py.ylabel('PSF Contrast')
return
[docs]def makecn2():
"""Reproduce the C_n^2 profile used in Guyon 2005 (Table 4).
:OUTPUTS:
z : 1D NumPy array
altitude of each layer, in meters
dz : scalar
the thickness of each layer, in meters.
cn2 : 1D NumPy array
C_n^2 (turbulence) profile of atmosphere.
:EXAMPLE:
import ao
z, dz, cn2 = ao.makecn2()
# Compute various moments
mu0 = (dz * cn2).sum()
mu2 = (dz * cn2 * z**2).sum()
mu53 = (dz * cn2 * z**(5/3.)).sum()
z2 = (mu2/mu0)**0.5
hbar = (mu53/mu0)**0.6
# Verify z_2 and hbar, in Table 4 of Guyon 2005:
print 'z_2 is: %1.1f km' % (z2/1000.)
print 'hbar is: %1.1f km' % (hbar/1000.)
"""
# 2012-12-05 20:11 IJMC: CReated
alts = np.array([.5, 1, 2, 4, 8, 16])*1e3
fracs = np.array([.2283, .0883, .0666, .1458, .3350, .1350])
dz = 100
z = np.arange(0,30e3,dz) # altitude, in m
cn2 = 0*z
for ii in range(alts.size):
ind = ((np.abs(z - alts[ii]))==(np.abs(z - alts[ii])).min()).nonzero()[0][0]
cn2[ind] = fracs[ii]
return z, dz, cn2
[docs]def calc_XY(z, dz, cn2, alpha_arcsec, lam_sci, lam_wfs=None):
"""Calculate X and Y of (Eq. 14 of Guyon 2005)
:INPUTS:
z : 1D NumPy array
Altitudes of given C_n^2 profile, in meters
dz : scalar or 1D NumPy array
Width of each layer of 'z', in meters; same size as 'z'
cn2 : 1D NumPy array
C_n^2 (turbulence) profile of atmosphere; same size as 'z'
alpha_arcsec : scalar or 1D NumPy array
Angular separation for calculation, in arcsec
lam_sci : scalar
Wavelength of science observations, in meters
lam_wfs : scalar
Wavelength of wavefront sensing, in meters (only needed for
calculating dX, dY terms).
:OUTPUTS:
X, Y (if lam_wfs is None)
X, Y, dX, dY (otherwise)
"""
# 2012-12-05 20:05 IJMC: Created
# Compute Spatial frequency; ensure it is an Array
f = alpha_arcsec / 206265. / lam_sci
if not hasattr(f, '__iter__'):
f = np.array([f])
if dz is None:
dz = np.diff(z)
dz = np.concatenate((dz[0], 0.5*(dz[1:] + dz[0:-1]), dz[-1]))
cn2dz = cn2*dz
denom = cn2dz.sum()
X = (cn2dz * (np.cos(np.pi*z*f.reshape(f.size, 1)**2 * lam_sci))**2).sum(1) / denom
Y = 1. - X
if lam_wfs is None:
ret = X, Y
else:
dX = (cn2 * dz * (np.cos(np.pi*z*f.reshape(f.size, 1)**2 * lam_sci) - \
np.cos(np.pi*z*f.reshape(f.size, 1)**2 * lam_wfs))**2).sum(1) / denom
dY = (cn2 * dz * (np.sin(np.pi*z*f.reshape(f.size, 1)**2 * lam_sci) - \
np.sin(np.pi*z*f.reshape(f.size, 1)**2 * lam_wfs))**2).sum(1) / denom
ret = X, Y, dX, dY
return ret
[docs]def contrast_uncorrected(lam_r0, r0, z, dz, cn2, alpha_arcsec, lam_sci, D, X=None, Y=None):
""" Compute PSF contrast from uncorrected atmospheric phase and
amplitude errors (C0 and C1)
:INPUTS:
lam_r0 : scalar
Wavelength at which 'r0' is specified, in meters
r0 : scalar
Fried parameter, in meters; specified at wavelength of 'lam_r0'
z : 1D NumPy array
Altitudes of given C_n^2 profile, in meters
dz : scalar or 1D NumPy array
Width of each layer of 'z', in meters; same size as 'z'
cn2 : 1D NumPy array
C_n^2 (turbulence) profile of atmosphere; same size as 'z'
alpha_arcsec : scalar or 1D NumPy array
Angular separation for calculation, in arcsec
lam_sci : scalar
Wavelength of science observations, in meters
D : scalar
Diameter of primary mirror, in meters.
X, Y: 1D NumPy Arrays
See :func:`calc_XY`
:OUTPUTS:
C0, C1: 2-tuple of 1D NumPy arrays
PSF contrast from uncorrected atmospheric phase and amplitude
errors (respectivey)
:NOTE:
Note that the coefficient "0.484" in Eqs. 16 and 17 of
Guyon 2005 has been corrected here to its intended value of
0.22*0.22=0.0484 (O. Guyon, private communication, Jan 2013).
"""
# 2012-12-05 19:56 IJMC: Created
# 2013-01-09 16:15 IJMC: Corrected C0/C1 coefficient, per O. Guyon
if X is None or Y is None:
X, Y = calc_XY(z, dz, cn2, alpha_arcsec, lam_sci)
alpha = alpha_arcsec / 206265. # Convert to radians
c01_term = 0.0484*np.pi**2 * lam_r0**2 * lam_sci**(5/3.) / \
(alpha**(11/3.) * D**2 * r0**(5/3.))
C0 = c01_term * X
C1 = c01_term * Y
return C0, C1
[docs]def contrast_residual(lam_r0, r0, z, dz, cn2, v, alpha_arcsec, F, beta_p, beta_a, lam_wfs, lam_sci, D, X=None, Y=None):
""" Compute PSF contrast from residual atmospheric phase and
amplitude errors (C2 and C3, of Guyon 2005)
:INPUTS:
lam_r0 : scalar
Wavelength at which 'r0' is specified, in meters
r0 : scalar
Fried parameter, in meters; specified at wavelength of 'lam_r0'
z : 1D NumPy array
Altitudes of given C_n^2 profile, in meters
dz : scalar or 1D NumPy array
Width of each layer of 'z', in meters; same size as 'z'
cn2 : 1D NumPy array
C_n^2 (turbulence) profile of atmosphere; same size as 'z'
v : scalar
Atmospheric windspeed, in meters per second.
alpha_arcsec : scalar or 1D NumPy array
Angular separation for calculation, in arcsec
F : scalar
Source brightness at 'lam_wfs', in photons per second per
meters squared.
beta_p : scalar or 1D NumPy array (same size as 'alpha_arcsec')
Wavefront sensor sensitivity to phase errors; equals 1.0 for
an ideal sensor.
beta_a : scalar or 1D NumPy array (same size as 'alpha_arcsec')
Wavefront sensor sensitivity to amplitude errors; equals 1.0
for an ideal sensor.
lam_wfs : scalar
Wavelength of wavefront sensor operation, in meters
lam_sci : scalar
Wavelength of science observations, in meters
D : scalar
Diameter of primary mirror, in meters.
X, Y: 1D NumPy Arrays
See :func:`calc_XY`
:OUTPUTS:
C2, C3 : 2-tuple of 1D NumPy arrays
PSF contrast from residual atmospheric phase and amplitude errors
"""
# 2012-12-05 19:56 IJMC: Created
if X is None or Y is None:
X, Y = calc_XY(z, dz, cn2, alpha_arcsec, lam_sci)
alpha = alpha_arcsec / 206265. # Convert to radians
numer = ((lam_wfs)**4 * (lam_r0 * v)**2)**(1./3.)
denom = (lam_sci**13 * (r0 * alpha)**5)**(1./9.) * (F**(1./3.) * D)**2
C2 = 0.7475 * (beta_p**4 * X)**(1./3.) * numer / denom
C3 = 0.7475 * (beta_a**4 * Y)**(1./3.) * numer / denom
return C2, C3
[docs]def contrast_chromaticity(lam_r0, r0, z, dz, cn2, alpha_arcsec, lam_wfs, lam_sci, D, X=None, Y=None, dX=None, dY=None, C0=None, C1=None, T=None, P=None, fco2=None, pph2o=None):
""" Compute PSF contrast from chromaticity of phase and amplitude
errors (C4 and C5)
:INPUTS:
lam_r0 : scalar
Wavelength at which 'r0' is specified, in meters
r0 : scalar
Fried parameter, in meters; specified at wavelength of 'lam_r0'
z : 1D NumPy array
Altitudes of given C_n^2 profile, in meters
dz : scalar or 1D NumPy array
Width of each layer of 'z', in meters; same size as 'z'
cn2 : 1D NumPy array
C_n^2 (turbulence) profile of atmosphere; same size as 'z'
alpha_arcsec : scalar or 1D NumPy array
Angular separation for calculation, in arcsec
lam_wfs : scalar
Wavelength of wavefront sensor operation, in meters
lam_sci : scalar
Wavelength of science observations, in meters
D : scalar
Diameter of primary mirror, in meters.
X, Y, dX, dY: 1D NumPy Arrays
See :func:`calc_XY`.
T, P, fco2, pph2o: scalars
Atmospheric parameters for calculating refractive index. If
set to None, the defaults in :func:`nsdata.nAir` will be used.
:OUTPUTS:
C4, C5, C6 : 3-tuple of 1D NumPy arrays
PSF contrast from chromaticity of: phase errors, amplitude
errors, and atmospheric refractivity.
"""
# 2012-12-05 19:56 IJMC: Created
if X is None or Y is None or dX is None or dY is None:
X, Y, dX, dY = calc_XY(z, dz, cn2, alpha_arcsec, lam_sci, lam_wfs=lam_wfs)
if C0 is None or C1 is None:
C0, C1 = contrast_uncorrected(lam_r0, r0, z, dz, cn2, alpha_arcsec, \
lam_sci, D, X=X, Y=Y)
C4 = C0 * dX / X
C5 = C1 * dY / Y
# Refractive indices:
n_wfs = nAir(lam_wfs*1e6, T=T, P=P, fco2=fco2, pph2o=pph2o)
n_sci = nAir(lam_sci*1e6, T=T, P=P, fco2=fco2, pph2o=pph2o)
C6 = C0 * ((n_sci - n_wfs) / (1. - n_sci))**2
return C4, C5, C6