% -*- mode: Noweb; noweb-code-mode: python-mode -*- \section{PXD file} \label{sec:pxd-file} <>= from numpy cimport ndarray from utilities cimport rtd, vector, cuIntArray, cuFloatArray, cuDoubleArray from rayTracing cimport coordinate_system, zernikeS, ZernikeS, Coordinate_system, Coordinates from source cimport Source, bundle, Bundle cdef extern from "gmtMirrors.h": <> <> @ \section{GMT mirrors} \label{sec:gmt-mirrors} \index{gmtMirrors!python!GmtMirrors} \subsection{PXD file} \label{sec:pxd-file-5} \subsubsection{Class definition} \label{sec:class-definition} <>= cdef class GmtMirrors: cdef: readonly float focal_plane_distance, focal_plane_radius public GMT_M1 M1 public GMT_M2 M2 public float pointing_error_zenith, pointing_error_azimuth public bint project_truss_onaxis public bint M2_baffle @ \subsection{PYX file} \label{sec:pyx-file} <>= cdef class GmtMirrors: """ A class container from GMT_M1 and GMT_M2 classes Parameters ---------- M1_radial_order : int, optionnal The largest radial order of the Zernike polynomials on M1 segments, default to 0 M2_radial_order : int, optionnal The largest radial order of the Zernike polynomials on M2 segments, default to 0 M1_mirror_modes : unicode, optional The modal basis on the M1 segments either ""zernike"" or ""bending modes"", default: ""zernike"" M1_N_MODE : int, optional The number of modes, default: 0 M2_mirror_modes : unicode, optional The modal basis on the M2 segments either ""zernike"" or "Karhunen-Loeve" (https://s3-us-west-1.amazonaws.com/gmto.rconan/KarhunenLoeveModes.bin), default: ""zernike"" M2_N_MODE : int, optional The number of modes, default: 0 Attributes ---------- M1 : GMT_M1 The GMT M1 CEO class M2 : GMT_M2 The GMT M2 CEO class sphere_radius : float The curvature radius of the ray tracing reference sphere See also -------- GMT_M1 : the class for GMT M1 model GMT_M2 : the class for GMT M2 model Source : a class for astronomical sources cuFloatArray : an interface class between GPU host and device data for floats Examples -------- >>> import ceo The mandatory parameters are the size of the pupil plane in meter or in pixel >>> gmt = ceo.GMT_MX(25.5,256) If more that one source (lets say 3) is going to be propagated through the telescope: >>> gmt = ceo.GMT_MX(25.5,256, N_SRC=3) A combination of Zernike polynomials can be applied to M1 and M2 segments by specifying the largest radial order on each mirror >>> gmt = ceo.GMT_MX(25.5,256, M1_radial_order=8, M2_radial_order=14) A source is propagated (geometrically) through the telescope with the following procedure: >>> src = ceo.Source("R",rays_box_size=25.5,rays_box_sampling=256,rays_origin=[0.0,0.0,25]) >>> gmt.propagate(src) and the wavefront phase is retrieved either as a 2D map cuFloatArray object with >>> gpu_ps2d = src.phase() or as a 1D vector with >>> gpu_ps1d = src.wavefront.phase() """ def __init__(self, int M1_radial_order=0, int M2_radial_order=0, unicode M1_mirror_modes=None, unicode M2_mirror_modes=None, int M1_N_MODE=0, int M2_N_MODE=0, dict M1_mirror_modes_data=None, dict M2_mirror_modes_data=None): self.M1 = GMT_M1(radial_order=M1_radial_order, mirror_modes=M1_mirror_modes, N_MODE=M1_N_MODE, mirror_modes_data=M1_mirror_modes_data) self.M2 = GMT_M2(radial_order=M2_radial_order, mirror_modes=M2_mirror_modes, N_MODE=M2_N_MODE, mirror_modes_data=M2_mirror_modes_data) self.focal_plane_distance = -5.830 self.focal_plane_radius = 2.197173 self.pointing_error_zenith = 0.0 self.pointing_error_azimuth = 0.0 self.project_truss_onaxis = False self.M2_baffle = False def __getitem__(self,key): if key=="M1": return self.M1 elif key=="M2": return self.M2 else: raise KeyError("Available keys are: M1 or M2") def __setitem__(self,key,value): if key=="M1": self.M1 = value elif key=="M2": self.M2 = value else: raise KeyError("Available keys are: M1 or M2") def propagate(self,Source src,str where_to="exit pupil", bint M2_is_a_stop=True): """ Propagate the Source object to the pupil plane conjugated to M1 Parameters ---------- src : Source The Source object where_to: char, optional Either "exit pupil" or "focal plane"; default: "exit pupil" M2_is_a_stop: boolean, optional Set if M2 acts also a stop; default: True """ cdef double [:] __zenith__, __azimuth__, theta_x, theta_y, offset_zenith, offset_azimuth __zenith__ = src._zenith_ __azimuth__= src._azimuth_ theta_x = src.theta_x - self.pointing_error_zenith*np.cos(self.pointing_error_azimuth) theta_y = src.theta_y - self.pointing_error_zenith*np.sin(self.pointing_error_azimuth) offset_zenith = np.hypot(theta_x,theta_y); offset_azimuth = np.arctan2(theta_y,theta_x); src.updateDirections(offset_zenith, offset_azimuth) src._c_source.reset_rays() #src.reset() if M2_is_a_stop: self.M2.blocking(src.rays) self.M1.trace(src.rays) if self.project_truss_onaxis: src.rays.gmt_truss_onaxis() if self.M2_baffle: src.rays.gmt_m2_baffle() self.M2.trace(src.rays) # src.sphere_distance # src.rays.to_sphere(self.sphere_radius,sphere_distance = src.sphere_distance) if where_to=="exit pupil": src.rays.to_sphere(focal_plane_distance=self.focal_plane_distance, focal_plane_radius=self.focal_plane_radius) src.updateDirections(__zenith__, __azimuth__) src.opd2phase() if where_to=="focal plane": src.rays.to_focal_plane(self.focal_plane_distance, self.focal_plane_radius) def reset(self, dict state=None): """ Reset M1 and M2 mirror segments to their original locations and shapes """ cdef str mirror self.pointing_error_zenith = 0.0 self.pointing_error_azimuth = 0.0 self.M1.motion_CS.reset() self.M1.modes.reset() self.M2.motion_CS.reset() self.M2.modes.reset() if state is not None: for mirror in state: self[mirror].__ixor__(state[mirror]) def __invert__(self): """ Reset the telescope to its original optical prescription """ self.reset() def update(self, dict state): """ Update M1 and M2 rigid body motions using dictionary input: {'M1':{'Rxyz':val1,...},'M2',{'Txyz':val2,...}} """ cdef str mirror for mirror in state: self[mirror].__ixor__(state[mirror]) def __ixor__(self, dict state): """ Update M1 and M2 rigid body motions using dictionary input: {'M1':{'Rxyz':val1,...},'M2',{'Txyz':val2,...}} """ cdef str mirror for mirror in state: self[mirror].__ixor__(state[mirror]) return self def dump_log(self, bytes filename): """ Writes the M1 and M2 segment commands log to a file Parameters ---------- filename : bytes The name of the file """ np.savez(filename, M1_rigid_body = self.M1.motion_CS.data_log, M1_zernike = self.M1.zernike.data_log, M2_rigid_body = self.M2.motion_CS.data_log, M2_zernike = self.M2.zernike.data_log) def segmentsWavefrontGradient(self, Source src): """ Computes the gradient of the phase of the wavefront on each segment Parameters ---------- src : Source The Source object Returns ------- cuFloatArray The wavefront gradient average in an array of size [N_SRCx14,1] """ cdef cuFloatArray sxy sxy = cuFloatArray( host_data=np.zeros((7*src.N_SRC*2,1) ) ) src._c_source.wavefront.segments_gradient_average(sxy._c_gpu.dev_data, sxy._c_gpu.dev_data+7, src.rays.L, src._c_source.rays.d__piston_mask) return sxy @property def state(self): cdef dict d d = {'M1': {'Rxyz':self.M1.motion_CS['Rxyz'].copy(), 'Txyz':self.M1.motion_CS['Txyz'].copy(), 'modes':np.copy(self.M1.modes.a)}, 'M2': {'Rxyz':self.M2.motion_CS['Rxyz'].copy(), 'Txyz':self.M2.motion_CS['Txyz'].copy(), 'modes':np.copy(self.M2.modes.a)}} return d @ \section{Bending modes} \label{sec:bending-modes} \index{gmtMirrors!python!bending\_modes} \subsection{PXD file} \label{sec:pxd-file-6} <>= cdef cppclass bending_modes: double *d__x_BM double *d__y_BM double *d__BM double *d__BMS int BM_N_SAMPLE double BM_radius void setup(int, int) void setupKL(int, int) void setupPolish(int, int) void cleanup() void load() void load_reg() void load_polish() void load_KL() void nearest_neighbor(rtd *, rtd *, rtd *, int , rtd , int) void bilinear(rtd *, rtd *, rtd *, int , rtd , int) void update(rtd *) @ \subsubsection{Class definition} \label{sec:class-definition-6} <>= cdef class BendingModes: cdef: bending_modes *_c_bending_modes public int n_mode, N_SURF, NPX public rtd[:,:] a @ \subsection{PYX file} \label{sec:pyx-file-6} <>= cdef class BendingModes: def __cinit__(self, int N_MODE=0, int N_SURF=1): self._c_bending_modes = new bending_modes() self._c_bending_modes.setup(N_MODE, N_SURF); self.n_mode = N_MODE self.N_SURF = N_SURF self.a = np.zeros((N_SURF, N_MODE), dtype=np.float64) def __dealloc__(self): self._c_bending_modes.cleanup(); def load(self,grid='irregular'): cdef: cuDoubleArray x_BM, y_BM, BM if grid=='regular': self._c_bending_modes.load_reg() BM = cuDoubleArray(shape=(self.n_mode, self._c_bending_modes.BM_N_SAMPLE*\ self._c_bending_modes.BM_N_SAMPLE)) BM._c_gpu.dev_data = self._c_bending_modes.d__BM self.NPX = 169 return BM else: self._c_bending_modes.load() x_BM = cuDoubleArray(shape=(1,self._c_bending_modes.BM_N_SAMPLE)) x_BM._c_gpu.dev_data = self._c_bending_modes.d__x_BM y_BM = cuDoubleArray(shape=(1,self._c_bending_modes.BM_N_SAMPLE)) y_BM._c_gpu.dev_data = self._c_bending_modes.d__y_BM BM = cuDoubleArray(shape=(1,self._c_bending_modes.BM_N_SAMPLE)) BM._c_gpu.dev_data = self._c_bending_modes.d__BM return (x_BM,y_BM,BM) def interpolate(self,int NI, double di, int k_mode=1,grid='irregular'): cdef: cuDoubleArray BMi cuDoubleArray partial_x_BMi cuDoubleArray partial_y_BMi BMi = cuDoubleArray(shape=(NI,NI), dev_malloc=True) partial_x_BMi = cuDoubleArray(shape=(NI,NI), dev_malloc=True) partial_y_BMi = cuDoubleArray(shape=(NI,NI), dev_malloc=True) if grid=='regular': self._c_bending_modes.bilinear(BMi._c_gpu.dev_data, partial_x_BMi._c_gpu.dev_data, partial_y_BMi._c_gpu.dev_data, NI, di, k_mode) else: self._c_bending_modes.nearest_neighbor(BMi._c_gpu.dev_data, partial_x_BMi._c_gpu.dev_data, partial_y_BMi._c_gpu.dev_data, NI, di, k_mode) return (BMi,partial_x_BMi,partial_y_BMi) <> def reset(self): """ Resets the bending modes coefficients to zero and update the bending modes surface """ self.a[:] = 0 self._c_bending_modes.update(&self.a[0,0]) @property def radius(self): return self._c_bending_modes.BM_radius @ with <>= def update(self, rtd[:,:] modes=None, **kwargs): """ Updates the bending modes surface based on the bending modes coefficients Parameters ---------- modes : ndarray The modal coefficients array in the $7\times [[N_MODE]]$ matrix form """ if modes is not None: self.a[:] = modes self._c_bending_modes.update(&self.a[0,0]) @ \section{Karhunen--Loeve modes} \label{sec:KL-modes} \index{gmtMirrors!python!KarhunenLoeve} \subsubsection{Class definition} \label{sec:class-definition-6} <>= cdef class KarhunenLoeve: cdef: bending_modes *_c_bending_modes public int n_mode, N_SURF public rtd[:,::1] a readonly cuDoubleArray S @ \subsection{PYX file} \label{sec:pyx-file-6} <>= cdef class KarhunenLoeve: def __cinit__(self, int N_MODE=0, int N_SURF=1): self._c_bending_modes = new bending_modes() self._c_bending_modes.setupKL(N_MODE, N_SURF); self.n_mode = N_MODE self.N_SURF = N_SURF self.a = np.zeros((N_SURF, N_MODE), dtype=np.float64) self.S = cuDoubleArray(shape=(self._c_bending_modes.BM_N_SAMPLE, self._c_bending_modes.BM_N_SAMPLE)) self.S._c_gpu.dev_data = self._c_bending_modes.d__BMS def __dealloc__(self): self._c_bending_modes.cleanup(); def load(self,grid='irregular'): cdef cuDoubleArray BM self._c_bending_modes.load_KL() BM = cuDoubleArray(shape=(1, self._c_bending_modes.BM_N_SAMPLE*\ self._c_bending_modes.BM_N_SAMPLE)) BM._c_gpu.dev_data = self._c_bending_modes.d__BM return BM def interpolate(self,int NI, double di, int k_mode=1,grid='irregular'): cdef: cuDoubleArray BMi cuDoubleArray partial_x_BMi cuDoubleArray partial_y_BMi BMi = cuDoubleArray(shape=(NI,NI), dev_malloc=True) partial_x_BMi = cuDoubleArray(shape=(NI,NI), dev_malloc=True) partial_y_BMi = cuDoubleArray(shape=(NI,NI), dev_malloc=True) self._c_bending_modes.bilinear(BMi._c_gpu.dev_data, partial_x_BMi._c_gpu.dev_data, partial_y_BMi._c_gpu.dev_data, NI, di, k_mode) return (BMi,partial_x_BMi,partial_y_BMi) <> def reset(self): """ Resets the bending modes coefficients to zero and update the bending modes surface """ self.a[:] = 0 self._c_bending_modes.update(&self.a[0,0]) @ \section{Polishing map} \label{sec:polish-map} \index{gmtMirrors!python!PolishMap} \subsubsection{Class definition} \label{sec:class-definition-6} <>= cdef class PolishMap: cdef: bending_modes *_c_bending_modes public int n_mode, N_SURF public rtd[:,::1] a readonly cuDoubleArray S @ \subsection{PYX file} \label{sec:pyx-file-6} <>= cdef class PolishMap: def __cinit__(self, int N_MODE=1, int N_SURF=1): self._c_bending_modes = new bending_modes() self._c_bending_modes.setupPolish(N_MODE, N_SURF); self.n_mode = N_MODE self.N_SURF = N_SURF self.a = np.zeros((N_SURF, N_MODE), dtype=np.float64) self.S = cuDoubleArray(shape=(self._c_bending_modes.BM_N_SAMPLE, self._c_bending_modes.BM_N_SAMPLE)) self.S._c_gpu.dev_data = self._c_bending_modes.d__BMS def __dealloc__(self): self._c_bending_modes.cleanup(); def load(self): cdef cuDoubleArray BM self._c_bending_modes.load_polish() BM = cuDoubleArray(shape=(1, self._c_bending_modes.BM_N_SAMPLE*\ self._c_bending_modes.BM_N_SAMPLE)) BM._c_gpu.dev_data = self._c_bending_modes.d__BM return BM def interpolate(self,int NI, double di, int k_mode=1): cdef: cuDoubleArray BMi cuDoubleArray partial_x_BMi cuDoubleArray partial_y_BMi BMi = cuDoubleArray(shape=(NI,NI), dev_malloc=True) partial_x_BMi = cuDoubleArray(shape=(NI,NI), dev_malloc=True) partial_y_BMi = cuDoubleArray(shape=(NI,NI), dev_malloc=True) self._c_bending_modes.bilinear(BMi._c_gpu.dev_data, partial_x_BMi._c_gpu.dev_data, partial_y_BMi._c_gpu.dev_data, NI, di, k_mode) return (BMi,partial_x_BMi,partial_y_BMi) <> def reset(self): """ Resets the bending modes coefficients to zero and update the bending modes surface """ self.a[:] = 0 self._c_bending_modes.update(&self.a[0,0]) @ \section{Modes} \label{sec:modes} \index{gmtMirrors!python!modes} \subsection{PXD file} \label{sec:pxd-file-7} <>= cdef cppclass modes: double *d__BM double *d__BMS int *d__s2b int N, BM_N_SAMPLE, N_SET, n_mode, N_MODE double BM_radius void setup(char *, int, int) void setup(int, double, int, int, int *, double *, int, int) void cleanup() void update(rtd *) void reset_modes(double *) @ \subsubsection{Class definition} \label{sec:class-definition-6} <>= cdef class Modes: cdef: modes *_c_modes public ndarray a rtd *a_ptr public cuDoubleArray M public cuDoubleArray S public cuIntArray s2b @ \subsection{PYX file} \label{sec:pyx-file-6} <>= cdef class Modes: """ A class representing modal basis Parameters ---------- filename : bytes The name of the file with the modal basis N : int The number of surfaces to generate from the modal basis n_mode : int The number of modes used to generate the surfaces Attributes ---------- a : rtd The basis coefficients as a [NXn_mode] array M : cuDoubleArray The modal basis array S : cuDoubleArray The surfaces array s2b : cuIntArray The surface to basis mapping indices See also -------- cuDoubleArray : an interface class between GPU host and device data for doubles cuIntArray : an interface class between GPU host and device data for ints """ def __cinit__(self, bytes filename=None, int N=7, int n_mode=0, int Ni=0, double L=0.0, int N_SET=0, int N_MODE=0, int[:] s2b=None , double[:] M=None): self._c_modes = new modes() if filename is None: self._c_modes.setup(Ni,L,N_SET,N_MODE,&s2b[0],&M[0],N,n_mode) else: self._c_modes.setup(filename,N,n_mode) self.a = np.zeros((self._c_modes.N, self._c_modes.n_mode), dtype=np.float64) if n_mode>0: self.a_ptr = self.a.data else: self.a_ptr = NULL self.M = cuDoubleArray(shape=(self._c_modes.BM_N_SAMPLE**2, self._c_modes.N_SET*self._c_modes.N_MODE), order='F') self.M._c_gpu.dev_data = self._c_modes.d__BM self.S = cuDoubleArray(shape=(self._c_modes.BM_N_SAMPLE**2,N), order='F') self.S._c_gpu.dev_data = self._c_modes.d__BMS self.s2b = cuIntArray(shape=N) self.s2b._c_gpu.dev_data = self._c_modes.d__s2b def __dealloc__(self): self._c_modes.cleanup(); def update(self, rtd[:,:] modes=None, **kwargs): """ Updates the bending modes surface based on the bending modes coefficients Parameters ---------- modes : ndarray The modal coefficients array in the $7\times [[N_MODE]]$ matrix form """ if modes is not None: self.a[:] = modes self._c_modes.update(self.a_ptr) def reset(self): """ Resets the bending modes coefficients to zero and update the bending modes surface """ self.a[:] = 0 self._c_modes.update(self.a_ptr) @property def radius(self): return self._c_modes.BM_radius @property def N_SAMPLE(self): return self._c_modes.BM_N_SAMPLE @property def N_SET(self): return self._c_modes.N_SET @property def N_MODE(self): return self._c_modes.N_MODE @property def n_mode(self): return self._c_modes.n_mode @n_mode.setter def n_mode(self,int value): self._c_modes.n_mode = value @property def modes(self): return self.M @modes.setter def modes(self,double[:] value): self._c_modes.reset_modes(&value[0]) @ \section{GMT M1} \label{sec:gmt-m1-1} \index{gmtMirrors!python!GMT\_M1} \subsection{PXD file} \label{sec:pxd-file-1} <>= cdef cppclass gmt_m1: double *d__x_BM double *d__y_BM double *d__BM int BM_N_SAMPLE void test_ray_tracing() <> void preset(bundle *, rtd ) @ where <>= int M_ID rtd D_assembly, D_full, D_clear, L, height, conic_c, conic_k float *d__segment_reflectivity rtd *d__conic_c rtd *d__conic_k coordinate_system aperture_CS coordinate_system conic_CS coordinate_system rigid_body_CS coordinate_system motion_CS coordinate_system TT_CS void setup() void setup( modes *) void setup( zernikeS *) void cleanup() void update(vector , vector ,int ) void reset() void trace(bundle *) void traceall(bundle *) void blocking(bundle *) void global_tiptilt(float , float ) void remove(int *, int ) void keep(int *, int ) void track(float *, float *, float *, int, int) void locate(float *, float *, float *, int, int) void update_conic_c(rtd *) void update_conic_k(rtd *) void set_reflectivity(float *) @ \subsubsection{Class definition} \label{sec:class-definition-1} <>= cdef class GMT_M1: cdef: gmt_m1 *_c_gmt_m12 <> @ @ with <>= cdef: readonly cuFloatArray segment_reflectivity public ZernikeS zernike public BendingModes BM public KarhunenLoeve KL public PolishMap PM public Modes B public Coordinate_system aperture_CS public Coordinate_system conic_CS public Coordinate_system rigid_body_CS public Coordinate_system motion_CS public Coordinate_system TT_CS public dict D public bytes mirror_modes_type public cuDoubleArray conic_curv public cuDoubleArray conic_const @ \subsection{PYX file} \label{sec:pyx-file-1} <>= cdef class DummyModes: def reset(self): pass # GMT M1 cdef class GMT_M1: <> def __cinit__(self, unicode mirror_modes=u"None", int radial_order=0, int N_MODE=0, dict mirror_modes_data=None, **kwargs): self._c_gmt_m12 = new gmt_m1() <> <> def test(self): self._c_gmt_m12.test_ray_tracing(); def preset(self, Bundle rays, float margin): self._c_gmt_m12.preset(rays._c_bundle, margin) def clocking(self,seg_id,theta): """ Decomposes mirror clocking into segment Txyx and Rxyz Parameters ---------- seg_id : any The segment number ID (1 base) theta : any The mirror clocking angle in radian Returns ------- (ndarray,ndarray) (Txyz,Rxyz) """ <> cdef double[:,:] rigid_body_CS_origin rigid_body_CS_origin = self.rigid_body_CS.origin[:] a = 0 b = 0 c = theta <> idx = seg_id-1 o = np.pi*(3-2*idx)/6 xt = self.L*np.cos(o+c) yt = self.L*np.sin(o+c) zt = self.rigid_body_CS.origin[idx,2] ut = np.array([xt,yt,zt])[:,None] <> return (origin,euler_angles) @ <>= def backward_transform(u,R,o): return np.dot(R,u) + o def forward_transform(u,R,o): return np.dot(R.T,u-o) @ <>= Ra = np.array( [[1, 0, 0,],[0,np.cos(a),-np.sin(a)],[0,np.sin(a),np.cos(a)]] ) Rb = np.array( [[np.cos(b), 0, np.sin(b)],[0,1,0],[-np.sin(b),0,np.cos(b)]] ) Rc = np.array( [[np.cos(c), -np.sin(c),0],[np.sin(c),np.cos(c),0],[0,0,1]] ) RTT = np.dot(Rc,np.dot(Rb,Ra)) @ <>= origin = np.ravel(forward_transform(ut,self.rigid_body_CS.R[idx], rigid_body_CS_origin[idx,:][:,None])) RM = np.dot(self.rigid_body_CS.R[idx].T,np.dot(RTT,self.rigid_body_CS.R[idx])) euler_angles = np.zeros(3) euler_angles[0] = np.arctan2(RM[2,1],RM[2,2]) euler_angles[1] = np.arcsin(-RM[2,0]) euler_angles[2] = np.arctan2(RM[1,0],RM[0,0]) @ <>= self.mirror_modes_type = mirror_modes.encode() if mirror_modes==b"zernike": self.zernike = ZernikeS(radial_order,N_SURF=7) self._c_gmt_m12.setup( self.zernike._c_zernikeS) elif mirror_modes==b"bending modes": self.BM = BendingModes(N_MODE,N_SURF=7) self._c_gmt_m12.setup( self.BM._c_bending_modes) elif mirror_modes==b"Karhunen-Loeve": self.KL = KarhunenLoeve(N_MODE,N_SURF=7) self._c_gmt_m12.setup( self.KL._c_bending_modes) elif mirror_modes==b"polishing map": print "!!polishing map!!" self.PM = PolishMap(1,N_SURF=7) self._c_gmt_m12.setup( self.PM._c_bending_modes) else: self._c_gmt_m12.setup() <>= self.conic_curv = cuDoubleArray(shape=(1,7),dev_malloc=False) self.conic_curv._c_gpu.dev_data = self._c_gmt_m12.d__conic_c self.conic_const = cuDoubleArray(shape=(1,7),dev_malloc=False) self.conic_const._c_gpu.dev_data = self._c_gmt_m12.d__conic_k self.mirror_modes_type = mirror_modes.encode('utf-8') if self.mirror_modes_type==b"None": self._c_gmt_m12.setup() elif self.mirror_modes_type==b"zernike": self.zernike = ZernikeS(radial_order,N_SURF=7) self._c_gmt_m12.setup( self.zernike._c_zernikeS) else: if mirror_modes_data is None: self.B = Modes(self.mirror_modes_type,7,N_MODE) else: self.B = Modes(N=7,n_mode=N_MODE,**mirror_modes_data) self._c_gmt_m12.setup( self.B._c_modes) self.segment_reflectivity = cuFloatArray(shape=(7,)) self.segment_reflectivity._c_gpu.dev_data = self._c_gmt_m12.d__segment_reflectivity self.aperture_CS = Coordinate_system(7) self.conic_CS = Coordinate_system(7) self.motion_CS = Coordinate_system(7) self.rigid_body_CS = Coordinate_system(7) self.TT_CS = Coordinate_system(1) self.aperture_CS.init( &(self._c_gmt_m12.aperture_CS) ) self.conic_CS.init( &(self._c_gmt_m12.conic_CS) ) self.motion_CS.init( &(self._c_gmt_m12.motion_CS) ) self.rigid_body_CS.init( &(self._c_gmt_m12.rigid_body_CS) ) self.TT_CS.init( &(self._c_gmt_m12.TT_CS) ) self.D = {} @ <>= def __dealloc__(self): self._c_gmt_m12.cleanup() def __ixor__(self, dict d): self.motion_CS.update(**d) self.modes.update(**d) return self def trace(self, Bundle rays, bint vignetting=True): """ Ray tracing to the mirror Parameters ---------- rays : Bundle A ray bundle used to propagate the light See also -------- Bundle: a class representing a bundle of light rays used for geometric propagation """ if vignetting: self._c_gmt_m12.trace(rays._c_bundle) else: self._c_gmt_m12.traceall(rays._c_bundle) def blocking(self, Bundle rays): """ Makes the mirror acting as a stop Parameters ---------- rays : Bundle A ray bundle used to propagate the light See also -------- Bundle: a class representing a bundle of light rays used for geometric propagation """ self._c_gmt_m12.blocking(rays._c_bundle) def update(self, list origin=[0.0,0.0,0.0], list euler_angles=[0.0,0.0,0.0], int idx=0): """ Updates the position of the mirror segments Parameters ---------- origin : list of float, optional The location of a segment center; default: [0,0,0] euler_angles : list of float, optional The tip, tilt and clock angles of a segment; default: [0,0,0] idx : int The segment index from 1 to 7 """ assert idx>0 and idx<8, "The segment index must be between 1 and 7!" self.motion_CS.origin[idx-1,:] = origin self.motion_CS.euler_angles[idx-1,:] = euler_angles self.motion_CS.update() def global_tiptilt(self, float tip, float tilt): """ Applies a global tip--tilt to M1 Parameters ---------- tip : float The x axis angle tilt : float The y axis angle """ self._c_gmt_m12.global_tiptilt(tip, tilt); def remove(self, list seg_ID): """ Remove some segments and keep the others Parameters ---------- seg_ID : int The ID numbers of the segment to remove """ cdef int[:] _seg_ID_ _seg_ID_ = np.array(seg_ID,dtype=np.int32,ndmin=1) self._c_gmt_m12.remove(&_seg_ID_[0], _seg_ID_.size) def keep(self, list seg_ID): """ Keep some segments and remove the others Parameters ---------- seg_ID : int The ID numbers of the segment to keep """ cdef int[:] _seg_ID_ _seg_ID_ = np.array(seg_ID,dtype=np.int32,ndmin=1) self._c_gmt_m12.keep(&_seg_ID_[0], _seg_ID_.size) def transform(self,list xyz, list abc, list segId=None): """ Decomposes mirror Txyz and Rxyz (origin is conic vertex) into segment Txyx and Rxyz Parameters ---------- seg_id : list The segment number ID (1 base) xyz : list The GCS translations abc : list The GCS rotations Returns ------- (ndarray,ndarray) (Txyz,Rxyz) """ <> cdef ndarray rigid_body_CS_origin rigid_body_CS_origin = self.rigid_body_CS.origin[:] rigid_body_CS_origin[:,2] = self.height - self.rigid_body_CS.origin[:,2] if segId is None: segId = range(7) else: segId = [x-1 for x in segId] a = abc[0] b = abc[1] c = abc[2] <> for idx in segId: ut = backward_transform(rigid_body_CS_origin[idx,:][:,None],RTT,np.transpose(np.array(xyz,ndmin=2))) <> self.motion_CS.origin[idx,:] = origin self.motion_CS.euler_angles[idx,:] = euler_angles self.motion_CS.update() def track(self, ndarray x, ndarray y, ndarray z, int segId): """ Coordinate transformation from motion CS to global CS Parameters ---------- x : ndarray The X axis coordinate y : ndarray The Y axis coordinate z : ndarray The Z axis coordinate segId : int The ID number of the segment """ cdef cuFloatArray gx, gy, gz cdef ndarray out = np.zeros((3,x.size)) gx = cuFloatArray( host_data = x ) gy = cuFloatArray( host_data = y ) gz = cuFloatArray( host_data = z ) self._c_gmt_m12.track(gx._c_gpu.dev_data, gy._c_gpu.dev_data, gz._c_gpu.dev_data, x.size, segId-1) out[0,:] = gx.host() out[1,:] = gy.host() out[2,:] = gz.host() return out def locate(self, ndarray x, ndarray y, ndarray z, int segId): """ Coordinate transformation from global CS to motion CS Parameters ---------- x : ndarray The X axis coordinate y : ndarray The Y axis coordinate z : ndarray The Z axis coordinate segId : int The ID number of the segment """ cdef cuFloatArray gx, gy, gz cdef ndarray out = np.zeros((3,x.size)) gx = cuFloatArray( host_data = x ) gy = cuFloatArray( host_data = y ) gz = cuFloatArray( host_data = z ) self._c_gmt_m12.locate(gx._c_gpu.dev_data, gy._c_gpu.dev_data, gz._c_gpu.dev_data, x.size, segId-1) out[0,:] = gx.host() out[1,:] = gy.host() out[2,:] = gz.host() return out property D_assembly: def __get__(self): return self._c_gmt_m12.D_assembly property D_full: def __get__(self): return self._c_gmt_m12.D_full property D_clear: def __get__(self): return self._c_gmt_m12.D_clear property L: def __get__(self): return self._c_gmt_m12.L property height: def __get__(self): return self._c_gmt_m12.height property conic_c: def __get__(self): return self._c_gmt_m12.conic_c property conic_k: def __get__(self): return self._c_gmt_m12.conic_k property modes: def __get__(self): """ if self.mirror_modes_type==u"zernike": return self.zernike elif self.mirror_modes_type==u"bending modes": return self.BM elif self.mirror_modes_type==u"Karhunen-Loeve": return self.KL elif self.mirror_modes_type==u"polishing map": return self.PM """ if self.mirror_modes_type is None: return DummyModes() elif self.mirror_modes_type==b"zernike": return self.zernike else: return self.B property M_ID: def __get__(self): return self._c_gmt_m12.M_ID property reflectivity: def __get__(self): return self.segment_reflectivity.host() def __set__(self, float[:] data): self._c_gmt_m12.set_reflectivity(&data[0]) @property def radius_error(self): return 1./self.conic_curv.host() - 1./self._c_gmt_m12.conic_c @radius_error.setter def radius_error(self, double[:] val): cdef int k cdef double[:] c = np.zeros(7) for k in range(7): c[k] = 1./(1./self._c_gmt_m12.conic_c + val[k]) self._c_gmt_m12.update_conic_c(&c[0]) @property def conic_error(self): return 1 - (1-self.conic_curv.host())/(1-self._c_gmt_m12.conic_c) @conic_error.setter def conic_error(self, double[:] val): cdef int k cdef double[:] kk = np.zeros(7) for k in range(7): kk[k] = 1 - (1 - val[k])*(1-self._c_gmt_m12.conic_k) self._c_gmt_m12.update_conic_k(&kk[0]) @ with <>= """ A class to represent GMT M1 or M2 segmented mirror Parameters ---------- radial_order : int, optional The radial order of the last Zernike polynomials, default to 0 mirror_modes : unicode, optional The modal basis on the segment either ""zernike"" or ""bending modes"" for M1 or "Karhunen-Loeve" (https://s3-us-west-1.amazonaws.com/gmto.rconan/KarhunenLoeveModes.bin) for M2, default: ""zernike"" N_MODE : int, optional The number of modes, default: 0 Attributes ---------- D_assembly : rtd The mirror assembly diameter D_full : rtd The segment full aperture diameter D_clear : rtd The segment clear aperture diameter L : rtd The distance from the optical axis to the center of the tilted peripheral segments conic_c : rtd The inverse of the mirror radius of curvature conic_k : rtd The conic parametert M_ID : int The mirror ID # zernike : ZernikeS The figure of the segments as Zernike surface object aperture_CS : Coordinate_system The segment aperture coordinate system conic_CS : Coordinate_system The segment conic coordinate system motion_CS : Coordinate_system The segment motion coordinate system rigid_body_CS : Coordinate_system The segment rigid body coordinate system TT_CS : Coordinate_system The mirror global tip-tilt coordinate system Examples -------- >>> import ceo >>> M1 = ceo.GMT_M1() With multiple sources: >>> M2 = ceo.GMT_M2() With Zernike modes >>> M1 = ceo.GMT_M1(25.5,101,radial_order=4) An M1 segment (#2) is displaced in x of 2 micron and tilted in y of 50mas with >>> import math >>> theta = 50e-3*math.pi/180/3600 >>> M1.update(origin=[2e-6,0,0],euler_angles=[0,theta,0],idx=2) A global tip-tilt of M1 is achieved with >>> M1.global_tiptilt(theta,-theta) M2 pointing neutral and coma neutral tip-tilt are realized with >>> M2.pointing_neutral(theta,-theta) >>> M2.coma_neutral(theta,-theta) """ @ \section{GMT M2} \label{sec:gmt-m2-1} \index{gmtMirrors!python!GMT\_M2} \subsection{PXD file} \label{sec:pxd-file-2} <>= cdef cppclass gmt_m2: <> void pointing_neutral(float, float) void coma_neutral(float, float) @ \subsubsection{Class definition} \label{sec:class-definition-5} <>= cdef class GMT_M2: cdef: gmt_m2 *_c_gmt_m12 <> @ \subsection{PYX file} \label{sec:pyx-file-2} <>= # GMT M2 cdef class GMT_M2: <> def __cinit__(self, unicode mirror_modes=u"None", int radial_order=0, int N_MODE=0, dict mirror_modes_data=None, **kwargs): self._c_gmt_m12 = new gmt_m2() <> <> def pointing_neutral(self, float tip, float tilt): """ Applies a pointing neutral global tip--tilt to M2 Parameters ---------- tip : float The x axis angle tilt : float The y axis angle """ self._c_gmt_m12.pointing_neutral(tip, tilt); def coma_neutral(self, float tip, float tilt): """ Applies a coma neutral global tip--tilt to M2 Parameters ---------- tip : float The x axis angle tilt : float The y axis angle """ self._c_gmt_m12.coma_neutral(tip, tilt); def clocking(self,seg_id,theta): """ Decomposes mirror clocking into segment Txyx and Rxyz Parameters ---------- seg_id : any The segment number ID (1 base) theta : any The mirror clocking angle in radian Returns ------- (ndarray,ndarray) (Txyz,Rxyz) """ <> cdef double[:,:] rigid_body_CS_origin rigid_body_CS_origin = self.rigid_body_CS.origin[:] a = 0 b = 0 c = theta <> idx = seg_id-1 o = np.pi*(3-2*idx)/6 + np.pi xt = self.L*np.cos(o+c) yt = self.L*np.sin(o+c) zt = self.rigid_body_CS.origin[idx,2] ut = np.array([xt,yt,zt])[:,None] <> return (origin,euler_angles) @ \section{Stereoscopic edge sensors} \label{sec:ster-edge-sens} \index{gmtMirrors!python!stereoscopic\_edge\_sensors} \subsection{PXD file} \label{sec:pxd-file-3} <>= cdef cppclass stereoscopic_edge_sensors: int N, N_DATA vector *v0 vector *dv0 vector *v vector *dv void setup(gmt_m1 *) void cleanup() void data() @ \subsubsection{Class definition} \label{sec:class-definition-2} <>= cdef class StereoscopicEdgeSensors: cdef: stereoscopic_edge_sensors *_c_stereoscopic_edge_sensors public Coordinates v0, v, dv0, dv @ \subsection{PYX file} \label{sec:pyx-file-3} <>= cdef class StereoscopicEdgeSensors: """ A class for the GMT stereoscopic edge sensors. The stereoscopic model measures the coordinates of the 3-dimension vector joining a pair of edge sensors. Parameters ---------- mirror : GMT_M1 The GMT M1 mirror Attributes ---------- v0 : Coordinates The coordinates of the edge sensors for a perfectly aligned telescope v : Coordinates The coordinates of the edge sensors for a disturbed telescope dv0 : Coordinates The coordinates of the vector joining a pair of edge sensors for a perfectly aligned telescope dv : Coordinates The coordinates of the vector joining a pair of edge sensors for a disturbed telescope See also -------- Coordinates : an interface between an array of device vectors and an array of host coordinates """ def __cinit__(self,GMT_M1 mirror): self._c_stereoscopic_edge_sensors = new stereoscopic_edge_sensors() self._c_stereoscopic_edge_sensors.setup(mirror._c_gmt_m12) self.v0 = Coordinates((self._c_stereoscopic_edge_sensors.N,3)) self.v0.v = self._c_stereoscopic_edge_sensors.v0 self.dv0 = Coordinates((self._c_stereoscopic_edge_sensors.N_DATA,3)) self.dv0.v = self._c_stereoscopic_edge_sensors.dv0 self.v = Coordinates((self._c_stereoscopic_edge_sensors.N,3)) self.v.v = self._c_stereoscopic_edge_sensors.v self.dv = Coordinates((self._c_stereoscopic_edge_sensors.N_DATA,3)) self.dv.v = self._c_stereoscopic_edge_sensors.dv def __dealloc__(self): self._c_stereoscopic_edge_sensors.cleanup() def data(self): """ Computes the edge sensor measurements """ self._c_stereoscopic_edge_sensors.data() @ \section{Lateral edge sensors} \label{sec:lateral-edge-sensor} \index{gmtMirrors!python!lateral\_edge\_sensors} \subsection{PXD file} \label{sec:pxd-file-4} <>= cdef cppclass lateral_edge_sensors: int N_DATA vector *A0 vector *A vector *B0 vector *B vector *k_cam vector *k_laser rtd *d__x rtd *d__y rtd *d__d void setup(gmt_m1 *) void setup(gmt_m1 *, rtd) void cleanup() void data() @ \subsubsection{Class definition} \label{sec:class-definition-3} <>= cdef class LateralEdgeSensors: cdef: lateral_edge_sensors *_c_lateral_edge_sensors public Coordinates A, A0, B, B0, k_cam, k_laser public cuDoubleArray x, y @ \subsection{PYX file} \label{sec:pyx-file-4} <>= cdef class LateralEdgeSensors: """ A class for the GMT lateral displacement edge sensors. It applies to the case of a laser aiming at a camera. The measurements are the laser spot location within the frame of the camera. Parameters ---------- mirror : GMT_M1 The GMT M1 mirror Attributes ---------- A : Coordinates The laser location coordinates in the motion coordinate system of the segment they are attached to. B : Coordinates The camera location coordinates in the rigid body coordinate system of the segment they are attached to. B0 : Coordinates The camera location coordinates in the global coordinate system of the segment they are attached to. k_cam : Coordinates The vector joining a pair of laser/camera in the rigid body coordinate system of the segment that the camera is attached to in the case of a perfectly aligned telescope. k_laser : Coordinates The vector joining a pair of laser/camera in the rigid body coordinate system of the segment that the camera is attached to in the case of a disturbed telescope. x : cuDoubleArray The x-axis coordinates of the laser spot y : cuDoubleArray The y-axis coordinates of the laser spot See also -------- Coordinates : an interface between an array of device vectors and an array of host coordinates cuDoubleArray : an interface class between GPU host and device data for doubles """ def __cinit__(self,*args,**kwargs): self._c_lateral_edge_sensors = new lateral_edge_sensors() def __init__(self,GMT_M1 mirror,*args,**kwargs): if len(args)>0: self._c_lateral_edge_sensors.setup(mirror._c_gmt_m12,args[0]) else: self._c_lateral_edge_sensors.setup(mirror._c_gmt_m12) self.A = Coordinates((self._c_lateral_edge_sensors.N_DATA,3)) self.A.v = self._c_lateral_edge_sensors.A self.A0 = Coordinates((self._c_lateral_edge_sensors.N_DATA,3)) self.A0.v = self._c_lateral_edge_sensors.A0 self.B0 = Coordinates((self._c_lateral_edge_sensors.N_DATA,3)) self.B0.v = self._c_lateral_edge_sensors.B0 self.B = Coordinates((self._c_lateral_edge_sensors.N_DATA,3)) self.B.v = self._c_lateral_edge_sensors.B self.k_cam = Coordinates((self._c_lateral_edge_sensors.N_DATA,3)) self.k_cam.v = self._c_lateral_edge_sensors.k_cam self.k_laser = Coordinates((self._c_lateral_edge_sensors.N_DATA,3)) self.k_laser.v = self._c_lateral_edge_sensors.k_laser self.x = cuDoubleArray(shape=(self._c_lateral_edge_sensors.N_DATA,1),dev_malloc=False) self.x._c_gpu.dev_data = self._c_lateral_edge_sensors.d__x self.y = cuDoubleArray(shape=(self._c_lateral_edge_sensors.N_DATA,1),dev_malloc=False) self.y._c_gpu.dev_data = self._c_lateral_edge_sensors.d__y def __dealloc__(self): self._c_lateral_edge_sensors.cleanup() def data(self): """ Computes the edge sensor measurements """ self._c_lateral_edge_sensors.data() @ \section{Distance edge sensors} \label{sec:dist-edge-sens} \index{gmtMirrors!python!distance\_edge\_sensors} \subsubsection{Class definition} \label{sec:class-definition-4} <>= cdef class DistanceEdgeSensors(LateralEdgeSensors): cdef: cuDoubleArray _d_ readonly double[:,::1] d0 public float error_1sig_per_m @ \subsection{PYX file} \label{sec:pyx-file-5} <>= import numpy as np cdef class DistanceEdgeSensors(LateralEdgeSensors): """ A class for the GMT distance edge sensors. It applies to the case of a laser aiming at a retroreflector. The measurements are the distances between the laser spot and the retroreflector. Parameters ---------- mirror : GMT_M1 The GMT M1 mirror height : rtd The laser and retroreflector are placed at a distance of +/- height from the z=0 plane of the segments Attributes ---------- d : cuDoubleArray The distance between the laser and the retroreflector See also -------- cuDoubleArray : an interface class between GPU host and device data for doubles """ def __init__(self,GMT_M1 mirror, rtd height=0.25, float error_1sig_per_m=0.0): super().__init__(mirror, height) self._d_ = cuDoubleArray(shape=(self._c_lateral_edge_sensors.N_DATA,1),dev_malloc=False) self._d_._c_gpu.dev_data = self._c_lateral_edge_sensors.d__d self.d0 = np.zeros((self._c_lateral_edge_sensors.N_DATA,1)) self.error_1sig_per_m = error_1sig_per_m def calibration(self): self.data() self.d0 = self._d_.host() property d: def __get__(self): #self.data() if self.error_1sig_per_m>0: return self._d_.host()-self.d0 + \ np.random.randn(self._c_lateral_edge_sensors.N_DATA,1) * \ self.error_1sig_per_m*self._d_.host() else: return self._d_.host()-self.d0 @ \subsection{TT7} \label{sec:tt7} \subsubsection{Diffractive} \label{sec:diffractive} \index{gmtMirrors!python!TT7} <>= from shackHartmann cimport ShackHartmann from atmosphere cimport AtmosphereAbstract cimport numpy as np cdef class TT7(ShackHartmann): cdef: list segID GmtMirrors _gmt_ public AtmosphereAbstract atm readonly double[:,:] _c0_, _c_ readonly np.ndarray frames <>= import numpy as np cdef class TT7(ShackHartmann): def __init__(self, GmtMirrors gmt, *args, **kwargs): ShackHartmann.__init__(self,1, *args, **kwargs) self._gmt_ = gmt self.atm = None self._c0_ = np.zeros((7,2)) self._c_ = np.zeros((7,2)) self.frames = np.zeros((self.N_PX_FRAME,self.N_PX_FRAME,7),dtype=np.float32) self.segID = range(1,8) def process(self): pass def get_measurement(self): return np.ravel(np.array(self._c_)) def get_measurement_size(self): return 14 def calibrate(self, Source gs): cdef int k cdef float[:] buf for k in self.segID: <> self._c_shackHartmann.analyze(gs._c_source) self.frames[:,:,k-1] = self.frame.host() buf = self.c.host().ravel() self._c0_[k-1,0] = buf[0] self._c0_[k-1,1] = buf[1] self._gmt_.M1.keep(self.segID) def propagate(self, Source gs): cdef int k, n7 cdef float[:] buf cdef float d7 if self.atm is None: for k in self.segID: <> self._c_shackHartmann.propagate(gs._c_source) self._c_shackHartmann.camera.readout(self.camera.exposureTime, self.camera.readOutNoiseRms, self.camera.nBackgroundPhoton, self.camera.noiseFactor) self.frames[:,:,k-1] = self.frame.host() self._c_shackHartmann.process() buf = self.c.host().ravel() self._c_[k-1,0] = buf[0] - self._c0_[k-1,0] self._c_[k-1,1] = buf[1] - self._c0_[k-1,1] else: n7 = gs.rays.N_L d7 = gs.rays.L/(n7-1) for k in self.segID: <> self.atm.ray_tracing(gs,d7,n7,d7,n7,gs.timeStamp) self._c_shackHartmann.propagate(gs._c_source) self._c_shackHartmann.camera.readout(self.camera.exposureTime, self.camera.readOutNoiseRms, self.camera.nBackgroundPhoton, self.camera.noiseFactor) self.frames[:,:,k-1] = self.frame.host() self._c_shackHartmann.process() buf = self.c.host().ravel() self._c_[k-1,0] = buf[0] - self._c0_[k-1,0] self._c_[k-1,1] = buf[1] - self._c0_[k-1,1] self._gmt_.M1.keep(self.segID) def analyze(self, Source gs, **kwargs): cdef int k, n7 cdef float[:] buf cdef float d7 if self.atm is None: for k in self.segID: <> self._c_shackHartmann.analyze(gs._c_source) self.frames[:,:,k-1] = self.frame.host() buf = self.c.host().ravel() self._c_[k-1,0] = buf[0] - self._c0_[k-1,0] self._c_[k-1,1] = buf[1] - self._c0_[k-1,1] else: n7 = gs.rays.N_L d7 = gs.rays.L/(n7-1) for k in self.segID: <> self.atm.ray_tracing(gs,d7,n7,d7,n7,gs.timeStamp) self._c_shackHartmann.analyze(gs._c_source) self.frames[:,:,k-1] = self.frame.host() buf = self.c.host().ravel() self._c_[k-1,0] = buf[0] - self._c0_[k-1,0] self._c_[k-1,1] = buf[1] - self._c0_[k-1,1] self._gmt_.M1.keep(self.segID) @ with <>= self.reset() self._gmt_.M1.keep([k]) gs.reset() self._gmt_.propagate(gs)