% -*- 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} <>= def BuildLinearActiveOptics(filename, gmt_prms, gs_tt7_prms, gs_wfs_prms, wfs_prms, includeBM=True): """ """ <> <> <> <> <> @ 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: <>= 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: <>= 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: <>= 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: <>= 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: <>= 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} <>= def Wavefronts(filename,science_src_prms, stroke=1e-6,includeBM=True): <> <> <> @ Lets first computes the wavefront (amplitude [[a0]] and phase [[ps0]]) for an ideally collimated telescope: <>= 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. <>= 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$. <>= 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. <>= 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]]. <>= class LinearActiveOptics(object): def __init__(self, model, wavefronts): """ """ <> <> <> <> <> def __del__(self): self.model.close() self.wavefronts.close() <> <> <> <> <> <> @ The function [[__wfe_rms__]] computes the wavefront error rms in the telescope 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: <>= 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}. <>= 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 <>= 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: <>= 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 <>= 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: <>= 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. <>= 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]]: <>= @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. <>= 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, <>= @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 <> @ \item the imagelet size, <>= @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 <> @ \item the read--out noise, <>= @property def wfs_ron(self): return self.wfs.ron @wfs_ron.setter def wfs_ron(self,value): self.wfs.ron = value <> @ \item the background noise, <>= @property def wfs_nPhBackground(self): return self.wfs.nPhBackground @wfs_nPhBackground.setter def wfs_nPhBackground(self,value): self.wfs.nPhBackground = value <> @ \item the temporal controller, a simple integrator is used here, it is specified with a dictionnary: [[{'T':exposure_time,'tau':latency,'g':gain}]] <>= @property def wfs_controller(self): return self.wfs.controller @wfs_controller.setter def wfs_controller(self,value): self.wfs.controller = value <> @ \end{itemize} The residual WFE RMS combining all the error sources is given by <>= @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: <>= 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} <>= 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} <>= 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 <> <> <> <> <>