% -*- mode: Noweb; noweb-code-mode: python-mode -*-
\section{Introduction}
\label{sec:introduction}
The active optics linear model is divided in three parts.
A model is build first using the [[BuildLinearActiveOptics]] function.
This function uses \emph{CEO} to compute the linear relationships between the degrees of freedom (DOF) of $M_1$ and $M_2$ segments and the sensors.
The wavefronts in the exit pupil of the telescope corresponding to the DOFs of the segments are computed with the function [[Wavefronts]] and saved separately.
Then, the [[LinearActiveOptics]] class uses the model and the exit pupil wavefronts to estimate the performance of the system as a function of a set of properties.
The [[LinearActiveOptics]] class is designed according to the reactive programming paragdime.
It means that an object of the class [[LinearActiveOptics]] will re--evaluate the whole system error budget as soon as a system property has been updated.
\section{[[BuildLinearActiveOptics]]}
\label{sec:buildl}
<<BuildLinearActiveOptics>>=
def BuildLinearActiveOptics(filename, gmt_prms,
gs_tt7_prms, gs_wfs_prms,
wfs_prms, includeBM=True):
"""
"""
<<CEO objects instanciation>>
<<SH-WFS interaction matrix>>
<<TT7 interaction matrix>>
<<TT7 system wavefront calibration>>
<<saving model>>
@
The input argments to the [[BuilLinearActiveOptics]] function are [[filename]]: a \emph{shelve} file where the model is saved and a set of dictionaries.
The dictionnaries contain parameters that correspond to the input arguments of the \emph{CEO} classes representing the components of the active optics model i.e.:
\begin{itemize}
\item the GMT, [[gmt_prms]] ([[ceo.GMT_MX]]),
\item the TT7 guide star, [[gs_tt7_prms]] ([[ceo.Source]]),
\item the SH--WFS guide stars, [[gs_wfs_prms]] ([[ceo.Source]]),
\item the SH--WFS, [[wfs_prms]] ([[ceo.GeometricShackHartmann]]).
\end{itemize}
The TT7 sensor is a unique sensor which parameters are hard--coded into \emph{CEO}.
The optional keywords [[includeBM]]=True$|$False specify if the bending modes are included in the AcO reconstructor.
The \emph{CEO} objects are instanciated first:
<<CEO objects instanciation>>=
gmt = ceo.GMT_MX(**gmt_prms)
wfs_gs = ceo.Source(**gs_wfs_prms)
wfs = ceo.GeometricShackHartmann(**wfs_prms)
tt7_gs = ceo.Source(**gs_tt7_prms)
tt7 = ceo.GeometricTT7()
@
\subsection{Interaction matrices}
\label{sec:interaction-matrices}
The interaction matrix between the SH--WFSs and the DOFs of $M_1$ and $M_2$ segments is computed next:
<<SH-WFS interaction matrix>>=
print("@(BuildLinearActiveOptics)> WFS CALIBRATION ...")
wfs_gs.reset()
gmt.reset()
gmt.propagate(wfs_gs)
wfs.calibrate(wfs_gs,0.)
C = gmt.AGWS_calibrate(wfs,wfs_gs,decoupled=True,
fluxThreshold=0.5,includeBM=includeBM,
filterMirrorRotation=True,
calibrationVaultKwargs={'nThreshold':[2]*6+[0],
'insertZeros':[None]*7})
@
The interaction matrix between the TT7 and the tip and tilt of $M_2$ segments is now computed:
<<TT7 interaction matrix>>=
print("@(BuildLinearActiveOptics)> TT7 CALIBRATION ...")
gmt.reset()
gmt.propagate(tt7_gs)
tt7.calibrate(tt7_gs)
Dtt7 = gmt.calibrate(tt7,tt7_gs,
mirror = 'M2',
mode='segment tip-tilt',
stroke=1e-6)
@
The interaction matrix between the TT7 and the DOFs of $M_1$ and $M_2$ segments is computed last:
<<TT7 system wavefront calibration>>=
print("@(BuildLinearActiveOptics)> TT7 calibration of observables ...")
stroke = [1e-6]*4
D = []
D.append( gmt.calibrate(tt7,tt7_gs,mirror='M1',mode='Rxyz',stroke=stroke[0]) )
D.append( gmt.calibrate(tt7,tt7_gs,mirror='M2',mode='Rxyz',stroke=stroke[1]) )
D.append( gmt.calibrate(tt7,tt7_gs,mirror='M1',mode='Txyz',stroke=stroke[2]) )
D.append( gmt.calibrate(tt7,tt7_gs,mirror='M2',mode='Txyz',stroke=stroke[3]) )
if includeBM:
D.append( gmt.calibrate(tt7,tt7_gs,mirror='M1',mode='bending modes',stroke=1e-6) )
if includeBM:
N_MODE = gmt_prms["M1_N_MODE"]
D_s = [ np.concatenate([D[0][:,k*3:k*3+3],
D[2][:,k*3:k*3+3],
D[1][:,k*3:k*3+3],
D[3][:,k*3:k*3+3],
D[4][:,k*N_MODE:(k+1)*N_MODE]],axis=1)
for k in range(7)]
else:
D_s = [ np.concatenate([D[0][:,k*3:k*3+3],
D[2][:,k*3:k*3+3],
D[1][:,k*3:k*3+3],
D[3][:,k*3:k*3+3]],axis=1)
for k in range(7)]
@
The data that composed a model are saved in a \emph{shelve} file object:
<<saving model>>=
print("@(BuildLinearActiveOptics)> Saving model to %s"%filename)
db = shelve.open(filename)
db['gmt'] = gmt_prms
db['gs_wfs'] = gs_wfs_prms
db['wfs'] = wfs_prms
db['gs_tt7'] = gs_tt7_prms
db['C'] = C
db['Dtt7'] = Dtt7
db['D_s'] = D_s
db.close()
@
\section{[[Wavefronts]]}
\label{sec:wavefronts}
The [[Wavefronts]] function computes, through ray tracing with \emph{CEO}, the wavefronts in the GMT exit pupil corresponding to all the rigid body motions of the $M_1$ and $M_2$ segments and to the 162 bending modes of each $M_1$ segment.
The input arguments to the functions are the name of the \emph{shelve} file [[filename]] where the wavefront are saved and the dictionary of parameters of the science object [[science_src_prms]] similar to the input arguments of the constructor of the [[ceo.Source]] class.
The optional keyword are
\begin{itemize}
\item [[includeBM]]=True$|$False specifying if the wavefronts associated with the bending modes need to be computed,
\item [[stroke]]$=10^{-6}$, the amplitude of each DOF.
\end{itemize}
<<Wavefronts>>=
def Wavefronts(filename,science_src_prms,
stroke=1e-6,includeBM=True):
<<collimated wavefront>>
<<ray tracing>>
<<wavefront ordering and saving>>
@
Lets first computes the wavefront (amplitude [[a0]] and phase [[ps0]]) for an ideally collimated telescope:
<<collimated wavefront>>=
gmt = ceo.GMT_MX(M1_N_MODE=162, M1_mirror_modes='bending modes')
src = ceo.Source(**science_src_prms)
src>>(gmt,)
gmt.reset()
+src
ps0 = src.wavefront.phase.host()
a0 = src.wavefront.amplitude.host()
idx = a0==1
pupil_mask = idx
piston_mask = src.rays.piston_mask[0]
@
For each DOF, a ray tracing from the source through the telescope and to the exit pupil is performed.
The collimated wavefront [[ps0]] is removed from the each DOF wavefront.
The residual wavefront within the the GMT pupil is extracted and saved as a vector in a column of the matrix [[O]] for the rigid body motion and of the matrix [[B]] for the bending modes.
<<ray tracing>>=
k = 0
N_MODE = gmt.M1.modes.n_mode
if includeBM:
m = 0
B = np.zeros((a0.sum(dtype=np.int32),7*N_MODE))
#B = np.zeros((sampling**2,7*N_MODE))
O = np.zeros((a0.sum(dtype=np.int32),84))
#O = np.zeros((sampling**2,84))
sys.stdout.write('@(Wavefronts)> ray tracing:\n')
for mirror in ['M1','M2']:
sys.stdout.write(' . %s : '%mirror)
for segId in range(7):
sys.stdout.write('.')
sys.stdout.flush()
for mode in ('Rxyz','Txyz'):
for axis in range(3):
gmt.reset()
_Q_ = np.zeros((7,3))
_Q_[segId,axis] = stroke
gmt[mirror].motion_CS.update(**{mode:_Q_})
+src
_a_ = a0*src.wavefront.amplitude.host()
_ps_ = _a_*(src.wavefront.phase.host() - ps0)
O[:,k] = _ps_[idx] /stroke
#O[:,k] = src.wavefront.phase.host()
k+=1
if mirror=='M1' and includeBM:
for l in range(N_MODE):
gmt.reset()
gmt.M1.modes.a[segId,l] = stroke
gmt.M1.modes.update()
+src
_a_ = a0*src.wavefront.amplitude.host()
_ps_ = _a_*(src.wavefront.phase.host() - ps0)
B[:,m] = _ps_[idx] /stroke
#B[:,m] = src.wavefront.phase.host()
m+=1
sys.stdout.write('\n')
sys.stdout.flush()
@
The wavefronts are ordered segment wise and for each segment in the following order: $[R_{1,xyz},T_{1,xyz},R_{2,xyz},T_{2,xyz},B_1]_k,\forall k=1,\dots,7$.
<<wavefront ordering and saving>>=
OO = {}
OO['M1'] = []
OO['M2'] = []
OO['BM'] = []
a = 6
b = 7*a
for segId in range(7):
OO['M1'] += [O[:,segId*a:a*(segId+1)]]
OO['M2'] += [O[:,b+segId*6:b+6*(segId+1)]]
if includeBM:
OO['BM'] += [B[:,segId*N_MODE:N_MODE*(segId+1)]]
if includeBM:
W = [np.concatenate((OO['M1'][k],OO['M2'][k],OO['BM'][k]),axis=1) for k in range(7)]
else:
W = [np.concatenate((OO['M1'][k],OO['M2'][k]),axis=1) for k in range(7)]
W[-1] = np.delete(W[-1],[2,8],axis=1)
print("@(Wavefronts)> Saving wavefronts to %s"%filename)
db = shelve.open(filename)
db['science_src'] = science_src_prms
db['pupil_mask'] = pupil_mask
db['piston_mask'] = piston_mask
db['N_MODE'] = N_MODE
db['W'] = W
db.close()
@
\section{[[ModelData]]}
\label{sec:modeldata}
The [[ModelData]] class encapsulates the dictionnary of parameters of the \emph{CEO} model components.
The constructor of the class takes as input argument the dictionnary.
The dictionary properties becomes the class attributes.
Keywords arguments in the class constructor are also class attributes.
<<ModelData>>=
class ModelData(object):
def __init__(self,data,**kwargs):
self._data_ = data
for key in kwargs:
self.__dict__[key] = kwargs[key]
def __getattr__(self,name):
if name in self._data_.keys():
return self._data_[name]
else:
return self.__dict__[name]
@
\section{[[LinearActiveOptics]]}
\label{sec:linearactiveoptics}
The input arguments of the [[LinearActiveOptics]] class are a model and a wavefront \emph{shelve} files created respectively with the functions [[BuildLinearActiveOptics]] and [[Wavefronts]].
<<LinearActiveOptics>>=
class LinearActiveOptics(object):
def __init__(self, model, wavefronts):
"""
"""
<<photometry>>
<<model loading>>
<<wavefronts loading>>
<<load interaction matrices>>
<<system matrices>>
def __del__(self):
self.model.close()
self.wavefronts.close()
<<WFE RMS exit pupil>>
<<WFS noise variance>>
<<SVD truncation threshold property>>
<<WFS mutable properties>>
<<number of photon>>
<<wavefront samples>>
@
The function [[__wfe_rms__]] computes the wavefront error rms in the telescope exit pupil:
<<WFE RMS exit pupil>>=
def __wfe_rms__(self,_C_):
X = self.W.dot(_C_.dot(self.W.T))
return np.sqrt(np.sum(X.diagonal())/self.W.shape[0])*1e9
@
The function [[__wfs_noise_variance__]] computes the variance of the WFS centroids due to photon, read--out and background noise:
<<WFS noise variance>>=
def __wfs_noise_variance__(self):
nPh = self.nPhoton(self.gs_wfs.photometric_band,self.gs_wfs.magnitude[0])
nPhLenslet = self.wfs.exposureTime*self.wfs.photoElectronGain*\
nPh*(self.wfs.d)**2
self.wfs.noise_variance = wfsNoise(nPhLenslet,self.wfs.spotFWHM_arcsec,
self.wfs.pixel_scale_arcsec,
self.wfs.N_PX_IMAGE/self.wfs.BIN_IMAGE,
nPhBackground=self.wfs.nPhBackground,
controller=self.wfs.controller)
@
\subsection{Photometry}
\label{sec:photometry}
The photometry is defined according to the reference document \cite{photo}.
<<photometry>>=
self.photometry = {"V": {"wavelength":0.550E-6,"zero_point": 8.97E9},
"R": {"wavelength":0.640E-6,"zero_point":10.87E9},
"I": {"wavelength":0.790E-6,"zero_point": 7.34E9},
"J": {"wavelength":1.125E-6,"zero_point": 5.16E9},
"H": {"wavelength":1.654E-6,"zero_point": 2.99E9},
"K": {"wavelength":2.179E-6,"zero_point": 1.90E9},
"Ks": {"wavelength":2.157E-6,"zero_point": 1.49E9},
"R+I":{"wavelength":0.715E-6,"zero_point":24.46E9}}
@ %def self.photometry
From the photometry data and the star magnitude, the number of photon $[m^{-2}.s^{-1}]$ is obtained with
<<number of photon>>=
def nPhoton(self,band,magnitude):
return self.photometry[band]["zero_point"]*pow(10.0,-0.4*magnitude)
@ %def nPhoton
\subsection{Interaction matrices}
\label{sec:interaction-matrices}
@
Data from the model and the wavefronts files are loaded into the class attributes:
<<model loading>>=
self.model = shelve.open(model)
pixel_scale_arcsec = \
self.photometry[self.model['gs_wfs']['photometric_band']]['wavelength']*\
self.model['wfs']['BIN_IMAGE']/self.model['wfs']['d']*180*3600/np.pi/\
self.model['wfs']['DFT_osf']
self.gmt = ModelData(self.model['gmt'],variate=0.0)
self.wfs = ModelData(self.model['wfs'],
pixel_scale_arcsec = pixel_scale_arcsec,
spotFWHM_arcsec = 2*pixel_scale_arcsec,
ron=0.0, nPhBackground=0.0,
controller=None,
noise_variance=0.0,
variate=0.0)
self.tt7 = ModelData(None,variate=0.0)
self.gs_wfs = ModelData(self.model['gs_wfs'])
@ %def self.model self.gmt self.wfs self.gs_wfs
and
<<wavefronts loading>>=
self.wavefronts = shelve.open(wavefronts)
self.N_MODE = self.model['gmt']['M1_N_MODE']
self.W = sprs.csr_matrix( np.hstack([X[:,:self.N_MODE-self.wavefronts['N_MODE']] \
for X in self.wavefronts['W']]) )
@ %def self.wavefronts self.W self.N_MODE
The interaction matrices are loaded from the model:
<<load interaction matrices>>=
self.C = self.model['C']
D_s = self.model['D_s']
#print D_s[0].shape
#print D_s[-1].shape
Dtt7 = self.model['Dtt7']
@ %def self.C
From the interaction matrices, the matrices that describes the linear behavior of the system are computed.
<<system matrices>>=
Mtt7 = LA.inv(Dtt7)
P2 = np.zeros((12+self.N_MODE,2))
P2[6,0] = 1
P2[7,1] = 1
X = [P2.dot(Mtt7[k*2:(k+1)*2,[k,k+7]]) for k in range(6)] \
+ [np.delete(P2,[2,8],axis=0).dot(Mtt7[-2:,[6,13]])]
self._s_Mtt7 = sprs.block_diag(X)
X = [D_s[k][[k,k+7],:] for k in range(7)]
self._s_Dtt7 = sprs.block_diag(X)
self._s_Qtt7 = sprs.eye(self._s_Mtt7.shape[0]) - self._s_Mtt7.dot(self._s_Dtt7)
self._s_Dwfs = sprs.block_diag(self.C.D)
self._s_L = sprs.eye(self._s_Dwfs.shape[1])*1e-9
self.threshold = 1e-9
self.NOISE_WFS = 0.0
self.noise_wfs_wfe_rms = 0.0
self.gmt.variate = np.random.randn(self._s_Q.shape[1],1)
self.tt7.variate = np.random.randn(self._s_Stt7.shape[1],1)
self.wfs.variate = np.random.randn(self._s_Swfs.shape[1],1)
#self.D_s[-1] = np.insert(self.D_s[-1],[2,7],0,axis=1)
@ %def self._s_Mtt7 self._s_Dtt7 self._s_Qtt7 self._s_Dwfs self._s_L self.threshold self.NOISE_WFS self.noise_wfs_wfe_rms self.gmt.variate self.wfs.variate
The matrix that depends on the threshold of the truncation of the pseudo--inverse of the interaction matrix of the WFSs are re--evaluated at each update of [[self.threshold]]:
<<SVD truncation threshold property>>=
@property
def threshold(self):
return self.C.threshold
@threshold.setter
def threshold(self,value):
self.C.threshold = value
print(self.C.nThreshold)
self._s_Mwfs = sprs.block_diag(self.C.M)
self._s_Swfs = self._s_Mwfs
self._s_Qwfs = sprs.eye(self._s_Mwfs.shape[0]) - self._s_Mwfs.dot(self._s_Dwfs)
self._s_Q = self._s_Qwfs.dot(self._s_Qtt7)
self._s_S = self._s_Qwfs.dot(self._s_Mtt7.dot(self._s_Dtt7)) + \
self._s_Mwfs.dot(self._s_Dwfs)
self._s_Stt7 = self._s_Qwfs.dot(self._s_Mtt7)
self.FITTING = self._s_Q.dot(self._s_L.dot(self._s_Q.T))
#self.fitting_wfe_rms = self.__wfe_rms__(self.FITTING)
@ %def self._s_Mwfs self._s_Swfs self._s_Qwfs self._s_Q self._s_S self._s_Stt7 self.FITTING self.fitting_wfe_rms
\subsection{WFS mutable properties}
\label{sec:wfs-mutable-prop}
The noise of the wavefront centroids depends of the GS magnitude, the spot size, the read--out noise, the background npoise and the property of the temporal controller.
All these properties can be set with the attributes of [[self.gs]] and [[self.wfs]].
However, setting them with the following attributes will trigger the computation of the WFS noise covariance matrix and the associated RMS WFE.
<<WFS noise error update>>=
self.__wfs_noise_variance__()
self.NOISE_WFS = self.wfs.noise_variance*(self._s_Swfs.dot(self._s_Swfs.T))
#self.noise_wfs_wfe_rms = self.__wfe_rms__(self.NOISE_WFS)
@
The mutable properties are
\begin{itemize}
\item the GS magnitude,
<<WFS mutable properties>>=
@property
def wfs_gs_mag(self):
return self.gs_wfs.magnitude[0]
@wfs_gs_mag.setter
def wfs_gs_mag(self,value):
self.gs_wfs.magnitude = [value]*3
<<WFS noise error update>>
@
\item the imagelet size,
<<WFS mutable properties>>=
@property
def wfs_spotFWHM_arcsec(self):
return self.wfs.spotFWHM_arcsec
@wfs_spotFWHM_arcsec.setter
def wfs_spotFWHM_arcsec(self,value):
self.wfs.spotFWHM_arcsec = value
<<WFS noise error update>>
@
\item the read--out noise,
<<WFS mutable properties>>=
@property
def wfs_ron(self):
return self.wfs.ron
@wfs_ron.setter
def wfs_ron(self,value):
self.wfs.ron = value
<<WFS noise error update>>
@
\item the background noise,
<<WFS mutable properties>>=
@property
def wfs_nPhBackground(self):
return self.wfs.nPhBackground
@wfs_nPhBackground.setter
def wfs_nPhBackground(self,value):
self.wfs.nPhBackground = value
<<WFS noise error update>>
@
\item the temporal controller, a simple integrator is used here, it is specified with a dictionnary: [[{'T':exposure_time,'tau':latency,'g':gain}]]
<<WFS mutable properties>>=
@property
def wfs_controller(self):
return self.wfs.controller
@wfs_controller.setter
def wfs_controller(self,value):
self.wfs.controller = value
<<WFS noise error update>>
@
\end{itemize}
The residual WFE RMS combining all the error sources is given by
<<WFS mutable properties>>=
@property
def residual_wfe_rms(self):
return self.__wfe_rms__(self.FITTING+self.NOISE_WFS)
@
\subsection{Wavefront samples}
\label{sec:wavefront}
Wavefront samples are evaluated with the function [[wavefrontSamples]] taking as argument the number of sample [[N_SAMPLE]], the segment errror RMS [[segment_rms]], the WFS error rms [[wfs_rms]] and the exponent of the units conversion factor:
<<wavefront samples>>=
def wavefrontSamples(self,N_SAMPLE=1,segment_rms=0,
wfs_rms=0,tt7_rms=0,units_exponent=0):
"""
"""
Y = np.zeros((self.W.shape[1],N_SAMPLE))
sampling = self.wavefronts['science_src']['rays_box_sampling']
pupil_mask = self.wavefronts['pupil_mask']
X = np.zeros((sampling**2,N_SAMPLE))
if self.gmt.variate.shape[1]!=N_SAMPLE:
self.gmt.variate = np.random.randn(self._s_Q.shape[1],
N_SAMPLE)
if self.wfs.variate.shape[1]!=N_SAMPLE:
self.wfs.variate = np.random.randn(self._s_Swfs.shape[1],
N_SAMPLE)
if self.tt7.variate.shape[1]!=N_SAMPLE:
self.tt7.variate = np.random.randn(self._s_Mtt7.shape[1],
N_SAMPLE)
Y = self._s_Q.dot(self.gmt.variate)*segment_rms
Y -= self._s_Mtt7.dot(self.tt7.variate)*tt7_rms
Y -= self._s_Swfs.dot(self.wfs.variate)*wfs_rms
Z = self.W.dot(Y)
X[pupil_mask.flatten(),:] = Z
X = X.reshape(sampling,sampling,N_SAMPLE)
u = 10**-units_exponent
return (u*np.squeeze(X),u*np.std(Z,axis=0))
@
\section{Wavefront sensor noise}
\label{sec:wavefr-sens-noise}
<<wavefront sensor noise>>=
def wfsNoise(nPhLenslet,spotFWHM_arcsec,pixelScale_arcsec,
nPxLenslet,ron=0.0,nPhBackground=0.0, controller=None):
def readOutNoise(ron,pxScale,nPh,Ns):
return (pxScale*ron/nPh)**2*Ns**4/12
closed_loop_noise_rejection_factor = 1.0
if controller is not None:
s = lambda nu : 2*1j*np.pi*nu
G = lambda nu, T, tau, g : -g*np.exp(s(nu)*(tau+0.5*T))*np.sinc(nu*T)/(s(nu)*T)
E = lambda nu, T, tau, g : 1.0/(1.0+G(nu,T,tau,g))
N = lambda nu, T, tau, g : g*np.exp(s(nu)*tau)/(s(nu)*T)
H = lambda nu, T, tau, g : N(nu,T,tau,g)/(1.0+G(nu,T,tau,g))
T = controller['T']
closed_loop_noise_rejection_factor = \
quad(lambda x: np.abs( H(x,**controller) )**2,0,0.5/T)[0]*T*2
print("@(wfsNoise)> Closed-loop noise rejection factor: %.4f"%\
closed_loop_noise_rejection_factor)
ron_var = ron**2 + nPhBackground
sigma_noise = closed_loop_noise_rejection_factor*\
( (1e3*spotFWHM_arcsec/np.sqrt(2*np.log(2)*nPhLenslet)*MAS2RAD)**2 + \
readOutNoise(np.sqrt(ron_var),
pixelScale_arcsec*ARCSEC2RAD,
nPhLenslet,nPxLenslet) )
return sigma_noise
def noiseRejectionFactor(controller):
s = lambda nu : 2*1j*np.pi*nu
G = lambda nu, T, tau, g : -g*np.exp(s(nu)*(tau+0.5*T))*np.sinc(nu*T)/(s(nu)*T)
E = lambda nu, T, tau, g : 1.0/(1.0+G(nu,T,tau,g))
N = lambda nu, T, tau, g : g*np.exp(s(nu)*tau)/(s(nu)*T)
H = lambda nu, T, tau, g : N(nu,T,tau,g)/(1.0+G(nu,T,tau,g))
T = controller['T']
return quad(lambda x: np.abs( H(x,**controller) )**2,0,0.5/T)[0]*T*2
@
\section{The python module}
\label{sec:main-class}
<<LinearActiveOptics.py>>=
import os
import shelve
import sys
import numpy as np
import numpy.linalg as LA
from scipy.integrate import quad
import scipy.sparse as sprs
import ceo
ARCSEC2RAD = np.pi/180/3600#ceo.constants.ARCSEC2RAD
MAS2RAD = ARCSEC2RAD*1e-3#ceo.constants.MAS2RAD
arcs2rad = ARCSEC2RAD
<<BuildLinearActiveOptics>>
<<Wavefronts>>
<<ModelData>>
<<LinearActiveOptics>>
<<wavefront sensor noise>>